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PREFACE 


J 

This  report  consists  of  three  volumes  which  present  the  theory  and 
application  of  a valuable  data  reduction  tool,  the  analysis  of  covariance. 
Volume  I introduces  the  analysis  of  covariance  as  a general  linear  model 
(GLM)  and  then  expands  the  model  to  incorporate  the  multivariate  case, 
unequal  sample  size,  and  missing  observations  on  the  response  variable. 
Volume  I also  covers  the  analysis  of  covariance  for  nonparametric  data. 

Volumes  II  and  III  were  prepared  by  the  Department  of  Statistics, 
Oklahoma  State  University,  Stillwater,  Oklahoma  74074,  under  Air  Force 
Contract  F08635-76-C-0154,  with  the  Air  Force  Armament  Laboratory, 

Armament  Development  and  Test  Center,  Eglin  Air  Force  Base,  Florida  32542. 
The  contract  dealt  with  the  development  and  programing  of  the  method- 
ology for  evaluating  multiple  variable  data  with  missing  observations 
on  dependent  and  independent  variables  by  the  analysis  of  covariance 
method.  The  methodology  also  covers  case  for  unequal  sample  size.  This 
work  was  begun  in  January  1976  and  completed  in  December  1976.  This 
is  Volume  III. 

This  report  has  been  reviewed  by  the  information  Officer  (01)  and 
is  releasable  to  the  National  Technical  Information  Service  (NTIS).  At 
NTIS  it  will  be  available  to  the  general  public,  including  foreign 
nations. 

This  technical  report  has  been  reviewed  and  is  approved  for 
publication. 
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INTRODUCTION 


Many  techniques  in  multivariate  analysis,  such  as  multiple 
regression,  principal  component  analysis,  and  canonical  correlation, 
assume  that  one  has  an  array  of  numbers  (Yj^,  Y^,  •••,  Y^  , Xi_,  X 

1, 


1 *2'  "p  M1  2 

, n on  which  various  operations  are  performed. 


-•  V1 

In  ordinary  multiple  linear  regression  with  a single  dependent  variable, 
p will  be  1;  in  a designed  experiment,  some  of  the  X's  will  be  treatment 
and  block  indicators  taking  only  the  values  0 or  1. 


In  various  experimental  environments,  it  is  not  unconmon  that  some 
of  the  Y values  and/or  X values  will  be  missing,  resulting  in  one  or  more 
incomplete  observations:  an  animal  dies,  a test  tube  breaks,  a piece  of 
equipment  malfunctions,  or  a measuring  instrument  fails  to  record  as 
intended . 

Most  computer  programs  currently  available  require  that  complete 
observations  only  be  supplied  for  analysis.  One  therefore  has  the  choice 
of  discarding  any  observations  on  which  at  least  one  of  the  variables  is 
missing,  or  using  an  estimation  technique  of  some  kind  to  fill  in  the 
missing  values  before  submitting  the  data  for  computer  analysis.  If 
missing  values  occur  frequently,  the  first  alternative  can  result  in 
a drastic  reduction  in  the  amount  of  data  available  for  analysis.  The 
second  alternative  presents  difficult  problems  in  estimation;  no 
technique  is  available  which  works  well  under  all  circumstances.  (See, 
for  example,  Haitovsky  (1968).) 

In  Section  II  the  results  of  a technique  developed  to  handle  the 
situation  are  discussed.  Missing  values  occur  among  either  dependent 
or  independent  variables,  but  not  in  the  block  or  treatment  design 
indicator  variables.  Full  details  of  the  theoretical  development,  an 
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extension  of  the  work  of  Kleinbaum  (1973),  are  given  in  Volume  I.  The  ' 

proposed  approach  requires  neither  the  discarding  of  incomplete  observa- 
tions nor  the  prior  estimation  of  missing  values. 

In  Section  III,  user  instructions  are  provided  for  a computer 
program  written  to  perform  an  analysis  of  a set  of  data  by  the  methods 
given  in  Section  II.  Section  IV  provides  a worked  example,  while  a 

flow  chart  and  program  listing  are  given  in  Appendices  A and  B,  ‘ 

respectively. 


SECTION  II 


REVIEW  OF  THE  THEORETICAL  FRAMEWORK 


The  theory,  for  a multivariate  analysis  of  covariance  model  in 
which  missing  observations  occur  among  the  dependent  and/or  independent 
variables , was  developed  in  Volume  I , resulting  in  a general  procedure 
which  yields  Best  Asymptotic  Normal  (BAN)  estimators  of  the  design  param- 
eters and  regression  slopes  and  which  allows  hypothesis  tests  based  on 
statistics  which  are  asymptotically  distributed  as  chi-square  variables. 
Equation  10,  described  in  Section  III  of  Volume  I,  is  used  for  the  compu- 
tation of  a consistent  and  unbiased  estimate  of  the  variance -covariance 

A 

matrix  E.  The  value  of  E thus  obtained  is  then  used  to  obtain  estimates 
of  the  design  parameters,  regression  slopes,  and  chi-square  statistics 
for  the  hypothesis  tests.  Although  the  formulae  used  to  obtain  these  BAN 
estimates  and  test  statistics  are  different  in  form  from  those  described 
in  Volume  I,  they  are  mathematically  equivalent  and  are  used  in  the  actual 
calculations  only  to  make  the  computations  manageable  on  an  electronic 
computer.  These  alternative  formulae  are  obtained  from  the  Matrix  Modi- 
fied representation  of  the  MGMAC  model  which  is  described  in  the  following 
paragraphs . 


To  obtain  the  general  form  of  the  Matrix  Modified  representation 
of  the  MGMAC  model,  assume  that  there  are  n experimental  units  and  p 


response  variates  V^, 
written  in  the  form 


V in  total  and  that  the  model  can  be 
P 


E(Y)  = Xa  + ZB 
Var  (Y)  = In  8 E 

where  Y is  an  nxp  matrix  composed  of  p-variate  responses  on  n individuals 
with  missing  values  recorded  as  blanks, 

X is  an  nxmx  known  design  matrix  of  rank  R(x)  = rx(<mx<n)  corres- 
ponding to  the  classificatory  variables  of  the  model, 


a is  an  m^xp  matrix  of  unknown  parameters 

Z is  an  nxmz  matrix  composed  of  concomitant  variables  with  missing 
values  recorded  as  blanks, 

B is  an  m^xp  matrix  of  unknown  concomitant  parameters. 

T,  = ((a  g))  is  a PXP  positive  definite  matrix  of  usually  unknown  parameters 
which  represents  the  variance-covariance  matrix  of  any  row  of  Y, 
and  A 8 B is  the  Kronecker  Product  of  the  matrices  A and  B. 

Then  replace  Z by  Z*  where  Z*  is  derived  from  Z by  augmenting  Z (with 

O1  s in  place  of  missing  values)  by  a matrix  G of  dimension  (nxt)  composed 

of  t columns  each  with  a one  in  the  row  position  corresponding  to  the  missing 

values  in  Z and  zeroes  elsewhere.  (Note:  Z has  dimension  (nxm*)  where 

z 

m*  = m^  + t) . Thus,  the  Matrix  Modified  representation  of  the  MGMAC  model 
may  be  written  in  the  form 

E(Y)  = Xa  + Z*(^)  (1) 

Var(Y)  = In  8 Z 

where  6 is  a (txp)  matrix  of  unknown  parameters  due  to  the  t missing  values 
in  Z.  Missing  values  in  Y are  then  handled  by  an  approach  suggested  by 
Srivastava  (3)  for  incomplete  multiresponse  designs  and  later  by 
Kleinbaum  (l)  for  the  MGLM  model.  For  simplification  of  notation, 
the  model  [Equation  (l)]  is  first  written  in  the  form 

E(Y)  = A*  £ 

Var(Y)  = In  BZ 

where  A* (nxm*)  = [X  | Z*] , m*  = mx+m*  and  £(m*xp)  = 

experimental  units  are  then  divided  into  u disjoint  sets  of  experimental 

units  S, , S_,  ...,  S with  n.  units  in  S..  On  each  unit  in  the  set  S. 

1 2 u j j j 

only  q . (<p)  responses  are  observed  (i.e.,the  remaining  p-q^  response  variates 
are  missing  in  S^).  Therefore,  the  MGMAC  model  is  given  by 


The  n 


i = 1,2, . . . , 


t(Y.)  = A*  £ B. 

3 3 3 

Var(Y.)  = In  8 B . ' Z B . , j = 1,2, ...,u 
J 3 j 

where  Y^(n^x  qj  is  the  matrix  of  observations  for  the  jbh  set 

A*(n.xm*)  is  the  design  matrix  for  the  set  S., 

3 J J 

B.(pxq.)  is  the  incidence  matrix  of  rank  R(B.)  = q.  for  the  j th 

3 3 3 3 

set  of  experimental  units.  It  consists  of  0's  and  l's  and  is  defined 

by  B.  = ((b...  ))  where 

J se 

b...  = I 1 if  variate  V is  the  rth  ordered  variate 

(.3  3 I s — 

j sr  I 

j measured  in  the  jth  experimental  set 
I 0 otherwise . 

It  is  also  assumed  that  Y and  Y ' are  independent  if  ' and 

the  rows  of  Y^  are  independent  and  distributed  as  a q^-variate  multinormal 

vector  with  variance-covariance  matrix  B.'  Z B.. 

3 3 

Based  on  the  above  Matrix  Modified  representation  of  the  MGMAC  model, 
the  following  results  for  estimation  and  hypothesis  testing  are  obtained. 
Theorem  1.  A BAN  estimator  for  an  estimable  linear  set  0=  H'C  is  given 


- - u „ - u - 1 

= H'i  = H'[  ^ 8^  (Bj  ' Z ) _1BJ  * 8A*'A*]  ^[Bj  (B^  ' ZB.)~l  H A*'] 


where  Z is  defined  in  Theorem  4 of  Volume  I and  ^ is  formed  by  stacking 
the  columns  of  Yj  underneath  each  other. 

Theorem  2.  The  Wald  Statistic  for  testing  the  hypothesis  Hq:  H’^  = 0 


is  given  by 


Wn  - (H'£)'  {H ' [ l BjCBj'EBj)"^'  8 AJ'Aj]"  H>  _1  . 

Note:  The  above  theorems  follow  easily  from  similar  theorems  by  Kleinbaum  (l) • 
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SECTION  III 

INSTRUCTIONS  FOR  USE  OF  THE  COMPUTER  PROGRAM 

In  this  section  are  outlined  the  steps  for  preparation  of  data  cards 
and  control  cards  necessary  for  execution  of  the  computer  program.  A 
worked  example  illustrating  these  ideas  and  the  resulting  computer  out- 
put is  presented  in  the  next  section. 


PREPARATION  OF  DATA  CARDS 

The  methodology  and  computer  program  presented  in  this  report  were 
designed  to  assist  in  analyzing  data  from  an  experiment  in  which  a com- 
plete observation  consists  of  a block  identifier,  a treatment  identifier, 
numerical  values  on  one  or  more  concomitant  variables  (covariables,  inde- 
pendent variables),  and  numerical  values  on  one  or  more  dependent  varia- 
bles (response  variables).  For  example,  a complete  observation  may  come 
from  a single  shot  down  a firing  range,  where  block  corresponds  to  metal 
type,  treatment  corresponds  to  projectile  shape,  the  covariables  are  such 
things  as  initial  projectile  weight  and  initial  projectile  velocity,  and 
the  dependent  variables  are  residual  projectile  weight,  residual  projectile 
velocity,  and  plug  weight. 

For  this  program  the  data  card(s)  corresponding  to  a complete 
observation  should  contain  a block  number,  a treatment  number,  the 
covariables,  and  the  dependent  variables,  as  well  as  any  identification 
the  user  desires.  The  order  in  which  these  occur  on  the  card  is  not 
important.  It  is  necessary,  however,  that  both  blocks  and  treatments  be 
numbered  sequentially  beginning  with  "l." 


The  data  card  corresponding  to 

ment  2,  may  look  like: 

the 

observation  on  block  3, 

treat- 

TESTj06_  SITE  A 

3 2 

113 

.8  680 

110.7  507 

^8.9^ 

Identification 

Blk  & Trt 

X's 

Y's 

'] 


( 


! 


C 


I 


I 

i 


t 


' 


The  user  is  free  to  choose  any  number  he  desires  to  correspond  to 
missing  values  among  the  Y's  or  X's.  This  number  will  then  be  punched 
in  the  data  card  for  any  X or  Y value  which,  for  reasons  reviewed  in 
Section  I,  was  not  recorded  when  the  experiment  was  conducted.  Of  course 
the  selected  number  must  not  be  the  same  as  any  of  the  X or  Y values 
occurring  in  the  data  set  being  analyzed  and  must  also  be  of  a magnitude 
which  allows  it  to  be  coded  in  the  number  of  columns  provided  for  the 
X's  and  Y's.  Suppose  that,  in  the  example  above,  both  113.8  and  507 
had  not  been  recorded  and  were  to  be  treated  as  missing  values.  If  the 
missing  value  code  selected  by  the  user  was  UU . U , the  card  would  have 
been  punched  as : 

TEST  0 6 SITE  A3  2 OMj.U  680  110.7  8.9 

CONTROL  CARDS 

V/e  now  present  in  sequence  a description  of  the  control  cards  which 
must  accompany  every  job. 

Control  Card  1 - Current  Data  Set  Information 


The  first  21  columns  of  this  card  consist  of  seven  three-digit 
(i.e.,  13)  fields  which  contain  in  order  the  values  of 


NP  - The  number  of  dependent  variables  (Y's)  NP  > 1 
NT  - The  number  of  treatments  NT  t 1 
NB  - The  number  of  blocks  NB  > 1 
NK  - The  number  of  covariates  (X's)  NK  > 1 
NN  - The  number  of  observations 


IDGT  - An  input  parameter  to  the  LPSDOR  subroutine  (which  computes  a 
generalized  inverse  of  a matrix).  The  elements  of  the  matrix 
are  assumed  to  be  correct  to  IDGT  places.  Since  thin  program 
computes  the  matrix  elements  in  double  precision,  we  suggest 
that  the  user  supply  IDGT=lU. 
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NM^SS  - The  number  of  missing  values  in  the  covariates  only.  Do  not 
include,  in  this  count,  r.  issing  values  in  the  dependent  vari- 
ables . 


The  next  8 columns  are  blank.  Beginning  in  column  30  punch 


Columns  30-39:  EPS  - A test  value  for  zero  which  is  used  in  the 

DMFGR  subroutine  which  calculates  the  rank  of 
a matrix.  Suggest  user  let  EPS  = 1.0E  - 07. 

Columns  1*0-149:  D - The  double  precision  missing  value  code  described 

in  an  earlier  subsection  of  this  section.  In 
that  subsection  we  used,  as  an  example,  a missing 
vulue  code  of  l*!*.!*.  In  this  case  we  would  code 
I)  = 111*. *4  DO. 


Columns  50-63:  These  columns  contain  seven  two-digit  (12)  fields 

indicating  the  desired  print  options  for  the  analysis. 
Each  12  field  contains  either  "00,”  which  requests 
that  the  printing  be  suppressed,  or  "01,"  which 
requests  that  the  printing  not  be  suppressed.  In 
order,  the  7 codes  refer  to 


Columns 

50-51 : 

Columns 

52-53: 

Columns 

5**-55 : 

Columns 

56-57 : 

Columns 

58-59: 

Columns 

60-6l : 

Columns 

62-63: 

Print  option  for  listing  of  MAC  model. 

Print  option  for  listing  of  GMAC . model. 

Print  option  for  listing  of  MGMAC  model. 

Print  option  for  listing  of  sigma  and  its  Inverse. 

Print  option  for  listing  of  the  matrix  modified 
model . 

Print  option  for  subgroups  of  observations  corres- 
ponding to  the  different  patterns  of  missing  values. 

Print  option  for  beta  values. 


With  regard  to  the  user's  choice  of  the  print  options,  we  suggest 
that,  for  most  purposes,  it  would  be  sufficient  to  code  columns  50  through 
63  with  "01000001000001."  This  provides  the  user  with  a listing  of  the 
MAC  model,  sigma  and  its  inverse,  and  the  beta  values  as  well  as  results 
of  hypothesis  tests. 

Control  Card  2 - Variable  Format  for  the  Input  Data  Set 


On  this  card  the  user  supplies  a FORTRAN  format  statement  which  speci- 
fies the  columns  in  which  the  block  number,  treatment  number,  dependent 
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variables  (Yis  i = 1,  NP),  and  covariables  (X  , j = 1,  •••,  NK) , in 

J 

that  order,  are  to  be  found.  The  block  and  treatment  numbers  must  be  read 
according  to  an  Iv  format,  while  the  dependent  and  independent  variables 
would  ordinarily  be  read  with  F w.d  or  D w.d  formats.  This  card  must 
begin  with  a left  parenthsis  (in  column  l)  and  end  with  a right  parenthesis 
and  contain  no  intervening  blanks.  For  our  example  in  the  first  subsection 
of  this  section,  the  variable  format  card  would  be 

(Tl8,  Il,6x,Il,Tl+2,F5.1,lX,F4.0,3X,Fl*.l,T30,F5.1,lX,Fl+.0). 

Control  Card  3 - Number  of  Hypothesis  Matrices  Being  Supplied 


This  is  the  simplest  of  all  the  control  cards  and  merely  states  how 
many  hypothesis  matrices  are  to  be  supplied  to  provide  basis  for  hypoth- 
esis tests  for  the  particular  job.  This  is  an  integer,  NUMHYP,  which  is 
punched,  right  justified,  in  the  first  5 columns  of  the  card.  Ordinarily 
this  number  will  be  less  than  10,  so  that  only  column  5 need  be  punched. 

Remaining  Control  Cards 


Next,  we  must  supply  one  group  of  cards  for. each  of  the  hypotheses 
to  be  tested.  The  first  card  provides  the  following  information  relative 
to  each  hypothesis: 


Columns  1-5:  The  number  of  rows,  NR,  in  the  hypothesis  matrix.  This 
is  an  integer,  right  justified  in  columns  1-5. 


Columns  6-25:  Any  alphameric  code  which  identifies  the  hypothesis 

being  tested.  This  is  simply  an  identification  which 
will  be  listed  with  the  output  in  the  hypothesis  test- 
ing section  of  the  printout. 

Columns  26-27:  "00"  if  the  user  does  not  desire  to  have  the  test 

statistic  evaluated  for  this  hypothesis  on  all  NP 
response  variables  simultaneously.  "01"  if  the  user 
does  desire  to  have  the  test  statistic  evaluated  for 
this  hypothesis  on  all  NP  response  variables  simul- 
taneously. 


Columns  28-29:  "00"  if  the  user  does  not  desire  to  have  the  test  sta- 

tistic evaluated  for  this  hypothesis  on  the  first 
response  variable  separately.  "01"  if  the  user  desires 
to  have  the  test  statistic  evaluated  for  this  hypothe- 
sis on  the  first  response  variable  separately. 


« 


i 


u 
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Columns  ((26+  "00"  is  the  user  does  not  desire  to  have  the  test 

2NP)-(27+2NP) ) : statistic  evaluated  for  this  hypothesis  on  the  NPth 
response  variable  separately.  "0] " if  the  user  does 
desire  to  have  the  test  statistic  evaluated  for  this 
hypothesis  on  the  NP^h  response  variable  separately. 

The  remaining  NR  cards  in  the  group  for  a particular  hypothesis  give 
the  coefficients  for  each  fov  of  the  hypothesis  matrix,  one  row  on  each  of 
IIP  cards.  These  coefficients  are  read  from  consecutive  four-digit  fields 
according  to  an  Fb.2  format.  The  number  of  coefficients  read  per  card  will 
be  equal  to  (NB+KT+NK),  with  the  exception  that  if  the  data  set  is  composed 
of  one  block  only  (or  one  treatment  only)  the  number  will  be  ( NB+NT+NK-1 ) . 
This  number  will,  of  course,  vary  from  one  problem  to  the  next.  The  number 
of  rows  in  the  hypothesis  matrix,  NR,  will  also  vary  from  one  hypothesis  to 
the  next  and  from  problem  to  problem. 


If  we  denote  by  NRi,  i = 1,  ••*,  NUMHYP,  the  number  of  rows  in  the  ith 

hypothesis  being  tested,  then  the  number  of  control  cards  after  control  card 
„ . , . NUMHYP 

3 will  be  i (NRi  + i). 

i=l 

The  method  of  constructing  the  hypothesis  matrices  will  be  demonstrated 
in  the  next  subsection. 

EXAMPLES  ILLUSTRATING  CONSTRUCTION  OF  HYPOTHESIS  MATRICES 


In  Section  II,  we  outlined  the  general  formulation  of  't.  <=  assumed  model 

with 

E(Y)  = Xct  + ZB 

as  the  underlying  mean  structure  for  the  data.  Here  a was  an  by  p matrix 
of  unknown  design  parameters,  and  B was  an  mz  by  p matrix  of  unknown  regres- 
sion coefficients  for  the  concomitant  variables.  Let  us  stack  the  a matrix 
on  the  B matrix  to  form  a matrix  B.  In  keeping  with  the  notation  used  in 
relation  to  our  computer  control  cards,  this  will  be  a matrix  with  row 
dimension  (NB+NT+NE)  and  column  dimension  NP. 

Denoting  the  block  parameters  by  b^j's,  the  treatment  parameters  by 
tjj's,  arid  the  regression  coefficients  by  p^j's,  we  could  think  of  the 
I!  matrix  us 
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bn 

b12 

\p 

bNB1 

bNB.,  ••• 

C11 

t12 

C1 

NP 

tNT1 

cnt2  *•' 

tNTNP 

Bn 

612 

B1 

. NP 

bnk1 

6nk2 

6nK»P 

This 

makes  it 

easy 

NB  rows  for  block  parameters 


NT  rows  for  treatment  paramters 


NK  rows  for  regression  coefficients 


matrices  to  test  hypotheses  of  the  form  H_:  HB  = 0W,  where  H is  a full 

U n 

rank  hypothesis  matrix  of  dimension  r(=  rank  of  H)  by  (NB  + NT  + NK)  and  0., 

M 

is  the  r by  NP  matrix  containing  all  zeroes.  For  example,  the  hypothesis 

of  no  treatment  differences  (i.e.,  t^  = = ...  = t^  for  all  j = 1,2, ...,NP) 

3 j j 

can  be  tested  with  an  H matrix  containing  (NT-1)  orthogonal  rows.  Similarly, 
the  hypothesis  of  no  influence  of  covariates  (S^  * 0,  i = 1,...,NK;  j = 1,...,NP) 
can  be  tested  with  a full  rank  hypothesis  matrix  of  row  dimension  NK.  Some 
specific  examples  will  now  be  given  to  illustrate  the  construction  of  H matrices. 
Example  1: 

Suppose  we  have  two  blocks,  two  treatments,  one  covariate,  and  three  response 
variables.  Here  NB  * NT  = 2,  NK  = 1,  and  NP  = 3.  Thus,  B can  be  written 

bll  b12  b13 

b21  b22  b23 

B ■ 

fcll  t12  fc13 

fc21  C22  t23 

_B11  812  B13  _ 

To  test  the  hypothesis  of  no  treatment  differences,  H will  be  the  1 by  5 


H - [0  0 1-1  0] 

11 


matrix 


Note  that  HB  = 0 is  equivalent  to 


[(tn  - 

hi>  (ti2  " 

t22)  (t13  " t23)l  = t°  0 °] • 

To  test  the  hypothesis  of  no  influence  of  the  covariate,  H will 

be  the  1 by  5 matrix 

1 

H = [0  0 0 

0 1]. 

i 

Note  that  HB  = 

0 is 

equivalent  to 

i 

f6ll 

612 

S„]  - [0 

0 0]. 

1 Example  2: 

Suppose  we  have 

two 

blocks,  four 

treatments,  two  covariateB,  and 

three 

response  variables. 

Here 

NB  = 2,  NT  = 

4,  NK  = 2,  and  NP  = 3.  Thus,  B 

can 

be  written 

! 

~bll 

bl2  b13 

b21 

b22  b23 

hi 

t12  *”  13 

B = 

hi 

C22  fc23 

hi 

**32  C33 

hi 

*42  h 3 

hi 

612  R13 

hi 

P22  Pj23 

To  test  the  hypothesis  of  no  overall  treatment  differences,  H will  be 

the  3 by  8 matrix 

* 

~0 

0 1 -1  C 

0 0 0 1 

1 

H = 

0 

0 11-2 

1 

0 0 0' 

♦ 

0 

0 1 1 1 -3  0 0 J. 

. 

To  test  the  hypothesis  of  no  overall  block  differences,  H will  be 

the 

: , 1 by  8 matrix 

H = 

[i  • 

-1  0 0 0 0 0 0]. 

12 

..  . 


79  01  22  019 


To  test  the  hypothesis  of  no  overall  effect  due  to  the  covariates, 
II  will  be  the  2 by  8 matrix 


H 


00000010 
0 0 0 0 0 0 0 


J. 


To  test  the  hypothesis  of  no  effect  due  to  the  first  covariate  only 
H will  be  the  1 by  8 matrix 

H = [0  0 0 0 0 0 1 0]. 

To  test  the  hypothesis  of  no  difference  between  treatments  one,  three 
and  four,  H will  be  the  2 by  8 matrix 

0 0 1 0 -1  0 0 0 

[00  10  1 -2  00 


H = 


Example  3. 

Suppose  we  have  three  different  treatments  in  the  experiment  but  only  one 
block  and  one  covariate.  Suppose,  in  addition,  that  there  are  two  response 

variables.  Here  NB  = 1,  NT  = 3,  NK  = 1,  and  NP  = 2.  However,  the  row 
dimension  of  the  B matrix  is  (NB  + NT  + NK  — 1)  or  4,  since  we  do  not  require 
any  block  parameters.  (Refer  to  the  discussion  in  the  Control  Card  sub- 
section above.)  B will  be  the  4 by  2 matrix 


tll 

C12 

t21 

*22 

tv 

31 

32 

Bll 

P12 

B = 


The  hypothesis  matrix  for  testing  "H^:  No  Treatment  Differences"  will 
be  the  2x4  matrix 


j 1 -1  0 0 

L1  1 "2  0. 
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The  complete  set  of  punched  cards  necessary  to  execute  a job  using 
the  source  program  with  a level  G FORTRAN  compiler  are  submitted  to  the 
card  reader  in  the  following  order: 


data  cam/; 


COftTl'M.  t'Mib 


f cghtkoi.  r: 


'hM)  J 


< 


L 


bl  • 


ehtpt  maiii 


//LKEC.SYSIH  CL  * 


AESEMELEf.  SOL'PCi  LtCK 


f//.EHSU 


YSIN  EC  * 


EXEC  ASKFC 

hOPTKAU  EGUPCt  CECK 


//  hiki.w.iii  dc  • 


f II  KXKC  Kf 'KrCfM.f; , |<KO  I ON . CO*'./OOK 
JOli  CAM/ 


The  form  of  the  initial  job  card  and  the  various  job  control  cards 
may  vary  from  one  installation  to  another . The  user  should  verify  the 
validity  of  these  with  systems  personnel  at  his  particular  computer 
center. 
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SECTION  IV 


A WORKED  EXAMPLE 


A sample  problem  will  now  be  presented  and  the  necessary  programming 
steps  illustrated.  This  data  set  was  constructed  to  represent  a set  of 
results  from  a series  of  test  shots  involving  two  different  metals  from 
which  projectiles  were  made  (the  blocks)  and  three  different  projectile 
shapes  (the  treatments).  Each  block/treatment  combination  was  replicated 


five  times,  resulting  in  a total  of  30  observations.  Two  concomitant 
variables  were  measured  on  each  projectile.  These  were  X^=Initial  Projectile 
Weight  and  X2=Initial  Projectile  Velocity.  The  three  dependent  variables 
recorded  were  = Residual  Projectile  Weight,  Y2  = Residual  Projectile 
Velocity  , and  Y^  = Plug  Weight. 

The  data  were  generated  using  the  model 


fijk  = bi  + "j  + 6lXlljk+  62X2ijk+  £ljk 


where  Y^k  represents  the  kth  replication  of  treatment  j in  block  i for 
any  one  of  the  three  response  variates.  For  our  problem  i = 1,2;  j = 1,2,3; 
k = 1,2,..., 5.  The  disturbances  were  chosen  from  a table  of  random  standard 
normal  deviates  (2).  The  values  used  for  the  model  parameters  were 
as  shown  in  the  following  B matrix. 


bl 

Y1 

0 

Y2 

0 

Y3 

0 

CM 

J* 

2 

20 

0 

C1 

10 

8 

0 

C2 

10 

8 

0 

C3 

20 

20 

4 

61 

8 

0 

.005 

P2 

0 

7 

.01 
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The  data  resulting  from  this  simulation  were: 


Observat ion 

Block 

Treatment 

X1 

X2 

Y1 

Y2 

Y3 

1 

1 

1 

110 

710 

98.0 

504.4 

7.3 

2 

1 

1 

110 

710 

97.5 

503.2 

5.2 

3 

1 

1 

112 

695 

98.2 

494.5 

8.3 

4 

1 

1 

111 

705 

100.0 

501.1 

8.2 

5 

1 

1 

112 

690 

100.1 

490.9 

8.1 

6 

1 

2 

110 

710 

98.4 

504.6 

8.3 

7 

1 

2 

109 

705 

95.5 

501.2 

8.8 

8 

1 

2 

107 

700 

94.6 

497.0 

10.5 

9 

1 

2 

111 

700 

99.3 

498.8 

9.5 

10 

1 

2 

112 

710 

98.9 

504.3 

7.5 

11 

1 

3 

116 

810 

110.7 

586.4 

8.1 

12 

1 

3 

115 

770 

112.4 

558.9 

8.5 

13 

1 

3 

116 

790 

113.3 

572.6 

8.9 

14 

1 

3 

116 

800 

110.6 

581.1 

8.7 

15 

1 

3 

117 

790 

115.3 

572.9 

7.4 

16 

2 

1 

109 

705 

98.3 

522.8 

6.6 

17 

2 

1 

112 

705 

101.0 

521.3 

8.3 

18 

2 

1 

111 

690 

99.1 

510.0 

6.3 

19 

2 

1 

110 

700 

99.8 

517.2 

8.3 

20 

2 

1 

111 

710 

103.1 

524.3 

8.4 

21 

2 

2 

112 

700 

100.9 

517.6 

5.8 

22 

2 

2 

112 

690 

100.5 

510.5 

7.2 

23 

2 

2 

114 

695 

103.1 

514.6 

6.2 

24 

2 

2 

112 

700 

99.5 

517.3 

7.4 

25 

2 

2 

113 

705 

102.0 

522.9 

8.8 

26 

2 

3 

115 

795 

114.2 

597.1 

12.3 

16 

27 

2 

3 

116 

800 

114.8 

598.6 

12.5 

28 

2 

3 

117 

795 

115.9 

598.0 

12.7 

29 

2 

3 

116 

790 

116.0 

593.0 

12.3 

30 

2 

3 

117 

805 

115.2 

602.4 

12.4 

For  the  purpose  of  a realistic  example  in  which  we  have  some 
missing  observations  among  both  the  dependent  and  independent  variables, 
we  consider  as  missing  the  values  of 

Y2  on  observation  1 

Y^  and  on  observation  U 

X2  on  observation  18 

Y2  and  Y^  on  observation  2h . 

PREPARATION  OF  DATA  CARDS 

First  we  must  select  a value  to  use  as  the  missing  value  code.  Any 
negative  number  would  suffice,  as  would  zero.  We  have  chosen  D-50.0D0,  a number 
intermediate  in  magnitude  between  Y^  and  the  other  X and  Y variables.  The 
data  cards  were  punched  as  follows: 


Columns  1-5 
Column  10 

Column  12 

Columns  13-16 
Columns  17-20 
Columns  22-25 
Columns  27-30 
Columns  31-31* 


The  identification  SAMPLE 
Block  number 
Treatment  number 
X^,  in  format  Fh.O 

Xg , in  format  FU.O 

Y^ , in  format  FU.l  (Decimal  not  punched) 

Yg,  in  format  Fl*.l  (Decimal  not  punched) 

Y^»  in  format  FU.l  (Decimal  not  punched) 


For  example,  the  first  data  card  was,  beginning  in  column  1: 
SAMPLE  1 1 110  710  0980  0500  073 
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I 

I' 

Note  that  the  value  of  598.0  has  been  replaced  by  50.0,  since  In  this 
example  we  are  treating  it  as  a missing  value. 

PREPARATION  OF  CONTROL  CARDS 

Control  Card  1: 

To  complete  this  card,  we  determine  that  for  the  present  problem 
we  have : , 

NP  =■=  3 

NT  *=  3 

NB  = 2 \ 

NK  = 2 

NN  = 30  ! 

NMISS  * 2 

To  provide  a thorough  look  at  the  printout  available  from  the  program,  we  ^ ! 

shall  take  every  print  option.  Using  the  suggested  values  for  EPS  and 
IDGT,  the  first  control  card  has  the  form,  beginning  in  column  1: 

003003002002030014002  1.00E-07  50.0D0  01010101010101 

< 

Control  Card  2 : 

For  this  card,  we  may  use  any  permissible  FORTRAN  format  statement  to 
read  in  the  block  number,  treatment  number,  Y's  and  X's.  One  possibility 
for  our  data  would  be 

(T10 , U,  IX,  II,  T21,  2F5.1,  F4.1,  T14,  F3.0,  IX,  F3.0) 

Control  Card  3 and  Subsequent  Control  Cards : 

For  this  problem  we  shall  test  four  hypotheses  to  illustrate  the  flexibility 

I 

of  this  aspect  of  the  program.' 

(1)  H : No  Block  Differences 
o 

With  only  two  blocks,  the  H matrix  consists  of  one  row  which  compares 
the  two  block  effects.  Recalling  that  for  this  problem  the  matrix  of  unknown 
parameters  is 


f 


our  H matrix  is  H * [l  -1  0 0 0 0 oj  . If  we  elect  the  overall 

test  as  well  as  the  test  on  response  variates  two  and  three  separately, 
the  necessary  control  cards  are: 

00001N0  BLK  DIFF  01000101 

0100-100 

(2)  H : No  Treatment  Differences 

o 

The  H matrix  for  constrasting  our  three  treatments  requires  two  rows: 

~0  0 .1  -1  0 0 o- 

H = 

_0  0 1 1 -2  0 0_ 

The  necessary  cards,  electing  the  overall  test  as  well  as  and  Y ^ 
separately,  are: 

00002N0  TRT  DIFF  01010100 

000000000100-100 

0000000001000100-200 

(3)  H^:  No  Difference  Between  Treatments  1 and  2 

If  we  elect  the  overall  test  as  well  as  the  test  on  each  variate 
separately,  the  necessary  cards  are: 

00001N0  DIFF  BETW  TR  1&2  01010101 
000000000100-100 


(4)  ll():  No  Effort  Due*  to  Covariates 

Electing  the  same  test  options  as  for  the  previous  example,  our  control 
cards  are 

00002N0  EFFECT  COVARS  01010101 
000000000000000000000100 
0000000000000000000000000100 

These  control  cards  for  hypotheses  mist  be  preceded  by  control  card 
3,  which  for  our  problem  will  have  a "4"  punched  in  column  5,  indicating 
that  four  hypotheses  are  to  be  tested. 

THE  PRINTED  OUTPUT  FROM  THE  PROGRAM 

The  program  output,  when  all  print  options  are  exercised,  has  essentially 
10  parts. 

(1)  Problem  Dimensions: 

This  is  the  first  page  of  the  printout.  It  provides  the  user  with  the 
current  program  dimensions  and  dimensions  for  his  particular  job,  as  well 
as  the  print  options  chosen  for  the  job. 

(2)  Listing  of  the  Input  Data: 

Each  observation  is  given  with  its  block  value,  treatment  value,  dependent 
variables  and  Independent  variables  in  that  order. 

(i)  The  MAC  Model  Listing: 

The  response  variables  and  complete  design  matrix  are  listed  for  the 
multivariate  analysis  of  covariance  model. 

(4)  The  GMAC  and  MGMAC  Model  Listing: 

These  are  given  for  each  response  variate  separately,  along  with  the  rank 
of  the  corresponding  design  matrix. 
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(5)  Sigma  and  Smoothed  Sigma 

The  sigma  matrix  submitted  to  the  subroutine  SMOOTH  is  printed 
along  with  the  smoothed  sigma  matrix  which  the  subroutine  returns. 

(6)  Sigma  and  Sigma  Inverse 

Sigma  and  its  inverse,  as  oonputed  by  the  LPSDOR  subroutine,  are  given. 

(7)  The  Matrix  Modified  Model 

This  is  listed,  showing  the  extra  columns  added  to  the  design  matrix 
to  account  for  the  missing  values. 

(8)  The  Groups  of  Observations  with  Different  Patterns  of  Missing  Values 
These  are  the  "S  " sets  referred  to  in  Section  II.  The  design  matrices 

associated  with  these  sets  are  used  in  the  compilation  of  the  beta  estimates  and 
the  Wald  Statistics. 

(9)  Beta  Values 

The  beta  matrix  of  dimension  (NB  + NT  + NK)  by  NP  is  given. 

(10)  Hypothesis  Test  Results 

The  hypothesis  matrix  and  associated  Wald  Statistic  are  given  for  each 
of  the  hypotheses  to  be  tested. 

The  listing  of  the  computer  printout  for  the  sample  problem  follows  in 
the  remaining  pages  of  this  section. 

It  should  be  pointed  out  that  it  is  normal  for  underflow  errors  to  occur, 
particularly  in  large  scale  analyses.  We  might  also  add  that  the  error  code 
from  LPSDOR  is  printed  each  time  we  call  that  subroutine;  a successful  call 
to  LPSDOR  results  in  an  error  code  of  zero. 
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listing  of  input  data 


1 1 

i 

98.00 

50.00 

7 .30 

110.00 

710.00 

) 1 

i 

97.50 

503.20 

5.20 

110. 00 

710.00 

i 

98.20 

494.50 

8.30 

112.00 

695.00 

1 

i 

100.00 

501.10 

8.20 

111.00 

705.00 

t 

i 

100.1  0 

490.90 

8.  10 

112.00 

690.00 

« 

2 

98.40 

504.60 

8.30 

110.00 

710.00 

1 

2 

95.50 

SO  1.20 

8.80 

109.00 

705.00 

I 

2 

94.60 

497.00 

10.50 

107.00 

700.00 

I 

2 

99.  30 

498.80 

9.  SO 

111.00 

700.00 

t 

2 

99.90 

50  4 .30 

7.50 

112.00 

710.00 

3 

110.70 

586.40 

o 

o 

• 

o 

116.00 

50.00 

1 

3 

112.40 

558.90 

8.50 

115.00 

770. 00 

1 

3 

11 3.30 

572.60 

8.90 

116.00 

790.00 

1 

3 

110.60 

581.  10 

8.70 

116.00 

eoo.oo 

, 1 
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111.00 

705.00 

a.  i o 

1.00 

0.0 

1.00 

0.0 

0.0 

112.00 

690.00 

8.30 

1.00 

0.0 

0.0 

1 .00 

0.0 

1 10.00 

710.00 

a.eo 

1.00 

0.0 

0.0 

1 .00 

0.0 

109.00 

705.00 

10.50 

1.00 

0.0 

0.0 

1.00 

0.0 

107.00 

700.00 

9.50 

1.00 

0.0 

0.0 

1.00 

0.0 

111.00 

700.00 

7.50 

1.00 

0.0 

0.0 

1.00 

0.0 

112.00 

710.00 

8.50 

1.00 

0.0 

0.0 

c.o 

1.00 

115.00 

770. 00 

a. 90 

1.00 

0.0 

0.0 

0.0 

1.00 

116.00 

790.00 

8.70 

1 .00 

0.0 

0.0 

0.0 

1.00 

116.00 

800.00 

7.40 

1.00 

0.0 

0.0 

0.0 

1.00 

1 17.00 

790.00 

6.60 

0.0 

1.00 

1.00 

0.0 

0.0 

109.00 

705.00 

8.30 

0.0 

1.00 

1.00 

0.0 

0.0 

112.00 

705.00 

6.30 

0.0 

1.00 

1.00 

0.0 

0.0 

111.00 

50.00 

8.30 

0.0 

1.00 

1 .00 

0.0 

0.0 

110.00 

700.00 

8.40 

0.0 

1.00 

1 .00 

0.0 

0.0 

111.00 

710.00 

5.80 

0.0 

1 .00 

0.0 

1.00 

0.0 

1 12.00 

700.00 

7.20 

0.0 

1.00 

0.0 

1 .00 

0.0 

112.00 

690.00 

6.20 

0.0 

1.00 

0.0 

1.00 

0.0 

114.00 

695.00 

8.80 

0.0 

1.00 

0.0 

1 .00 

0.0 

113.00 

705.00 

12.30 

0.0 

1.00 

0.0 

0.0 

1.00 

115.00 

795.00 

12.50 

0.0 

1.00 

0.0 

0.0 

1.00 

116.00 

800.00 

12.70 

0.0 

1.00 

0.0 

0.0 

1.00 

117.00 

795.00 

12.30 

0.0 

1.00 

0.0 

0.0 

1.00 

116.00  . 

790.00 

12.40 

0.0 

1.00 

0.0 

0.0 

i.oo 

117.00 

805.00 

THE  VALUES  OF  V ANO  A FOR  VARIATE  3 FOLLC* 


MCMAC  MODEL 


7.30 

1.00 

0.0 

1.00 

0.0 

0.0 

110.00 

710.00 

0.0 

5.20 

1.00 

0*0 

1 .00 

0.0 

0.0 

110.00 

710.00 

0.0 

8.30 

1.00 

0.0 

1.00 

0.0 

0.0 

112.00 

695.00 

0.0 

8 « 20 

1 .00 

0.0 

1 .00 

0.0 

0.0 

111.00 

705.00 

0.0 

8.10 

1.00 

0.0 

1.00 

0.0 

0.0 

112.00 

690.00 

0.0 

8.30 

1.00 

0.0 

0.0 

1.00 

0.  0 

110.00 

710.00 

0.0 

8.80 

1.00 

0.0 

0.0 

1 .00 

0.0 

109.00 

705.00 

0.0 

10.50 

1 .00 

0.0 

0.0 

1 .00 

0.0 

107.00 

700.00 

0.0 

9.50 

1 .00 

0.0 

0 .0 

1.00 

0.0 

111.00 

700.00 

0.0 

7.50 

1.00 

0.0 

0.0 

1 .00 

0.0 

112.00 

710.00 

0.0 

8.50 

1.00 

0.0 

0.0 

0.0 

1. 00 

1 IS. 00 

770.00 

0.0 

8.90 

1 .00 

0.0 

0.0 

0.0 

1.00 

1 16.00 

790.00 

0.0 

8.70 

1.00 

0.0 

0.0 

0.0 

1.00 

116.00 

800.00 

0.0 

7.40 

1.00 

0.0 

0 .0 

0.0 

1.00 

117.00 

790.00 

0.0 

6.60 

0.0 

1.00 

1 .00 

0.0 

0.0 

109.00 

705.00 

0.0 

8.30 

0.0 

1.00 

1 .00 

c.o 

0.0 

112.00 

705.00 

0.0 

6.30 

0.0 

1 .00 

1 .00 

0.0 

0.0 

111.00 

0.0 

1.00 

8.30 

0.0 

1.00 

1.00 

0.0 

0.  0 

110.00 

700.00 

0.0 

8.40 

0.0 

1.00 

1 .00 

0.0 

0.0 

111.00 

710.00 

0.0 

5.80 

0.0 

1.00 

0.0 

1 .00 

0.0 

112.00 

700.00 

0.0 

7.20 

0.0 

1.00 

0.0 

1.00 

0.0 

112.00 

690.00 

0.0 

6.  20 

0.0 

1.00 

0.0 

1.00 

0.0 

114.00 

695.00 

0.0 

8.80 

0.0 

&.00 

0.0 

1.00 

0.0 

113.00 

705.00 

0.0 

12.30 

0.0 

1.00 

0.0 

0.0 

1.00 

115.00 

795.00 

0.0 

12.50 

0.0 

1.00 

0.0 

0.0 

1.00 

1 16.00 

800.00 

0.0 

12.70 

0.0 

1 .00 

0.0 

0.0 

1.00 

117.00 

795 . 00 

0.0 

12.30 

0.0 

1.00 

0.0 

0.0 

1.00 

116.00  . 

790.00 

0.0 

12.40 

0.0 

1.00 

0.0 

0.0 

1.00 

117.00 

805.00 

0.0 

RANK  OF 

THE  DESIGN 

MATRIX 

FOR  VARIATE  3 

FOR  THE 

MGMAC  MODEL  IS  5 

. 

: 


ERROR  CODE  FROM  LPSDCR 


0 


1 


INPUT  MATRIX  TO 
1.0201164 
-.11701501 
.19631602 


SMOOTH 

-.11701501 
.69946443 
. 1 1176522 


• 19631602 
.11176522 
2.2082109 


MATRIX  OUTPUT  QY 
1.0281164 
-.1170 1501 
. 19631602 


SMOOTH 

-.11 701 501 
.69946443 
.1 l 176522 


. 19631602 
• 111  76522 
2.2882109 


THE  VALUE  OF  SIGMA  FOLLOWS 

1.02812  -.117015 

-.117015  .699464 

•196316  .111765 

ERROR  CODE  FROM  LPSDOR 

I 

the  value  of  sigma  inverse  follows 


• 196316 

• 1 1 1765 
2.28821 


i 


l 

' 


1.01196  .184607 

.184607  1 .47459 

-.9583790-01  -.8786300-01 


-. 9583790-01 
-.8786300-01 
• 449537 
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700.00 


••••  THERE  ARE  7 DIFFERENT  POSSIBLE  GROUPS  KUMBEaFD  FROM  l TO  7 THOUGH  IN  GENERAL  NOT  ALL  GROUPS  WILL  APPEAR 


ERROR  COOE  FROM 


LPSOOB  0554  625EAFS6  005FIA20  00000008  005912C4 


**•**•  HYPOTHESIS  TESTING  SECTION 


* HYPOTHESIS  MATRICES  SHALL  BE  USED. IN  TURN, FOR  COMPUTING  CHI  SQUARE  STATISTICS 


EACH  MATRIX  SHOULD  HAVE 


7 COLUMNS  AS  FOLLOWS  J 


THE  FIRST 
THE  NEXT 

the  last 


2 COLUMNS  CORRESPOND  TO  BLOCK  PARAMETERS 

3 COLUMNS  CORRESPOND  TO  TREATMENT  PARAMETERS 

2 COLUMNS  CORRESPOND  TO  CCVARIATE  COEFFICIENTS 


• ••«  LISTING  OF  HYPOTHESIS  MATRIX  I FOLLOWED  BY  ITS  EXTENSION  FOR  MODIFIED  MODEL. HYPOTH  ID  IS:  NO  BLK  D1FF 


• •••  listing  a r HYPOTHESIS  MATRIX  2 FCLLCxED  by  its  extension  for  modified  model. mvpoth  id  is:  no  trt  oiff 


••*•••  ***«.ts  op  hypothesis  rests  on  ikoivioual  variates  for  hvpoth  MATRIX  2 WITH  to:  NO  TRT  OIFF 


••••  LISTING  of  hypothesis  matrix  a fclloweo  by  its  extension  for  mooifieo  model. hvpotm  10  is:  no  effect  covars 


RESULTS  OF  HYPOTHESIS  TESTS  ON  IKOIVIOUAL  VARIATES  FOR  HVROTH  MATRIX  4 WITH  IOS  NO  EFFECT  COVARS  ••*••• 


«REES  OF  FREEDOM 
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3.  Haitovsky,  Y.  (1968).  Missing  Data  in  Regression  Analysis,  J.  Roy. 
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APPENDIX  A 
ITXJW  Q1AKT 


READ  INPUT 
PARAMETERS 


READ  INPUT  DATA 
BLOCK  & TRT  CODES,  DEPENDENT y 
AND  INDEPENDENT  VARIABLES 


l 


WRITE  INPUT  DATA 


7 


SET  UP  BLOCK  AND  TRT 
DESIGN 


/ BUILD  COVARIATE 
\ MATRIX 


BUILD  DESIGN  MATRIX 
AND  AUGMENT  WITH  COVARIATE 
MATRIX 


\ 


\ 


WRITE  THE  MATRIX 
FORM  OF  THE  MAC  MODEL 


' CONSTRUCT  THE  VARIATE  FORM  OF  THE  MODEL 
REMOVING  MISSING  VALUES  FROM  THE  DEPENDENT 
VARIABLES  AND  THE  CORRESPONDING  ROWS  OF 
THE  DESIGN  MATRIX 


WRITE  THE  VARIATE  FORM  OF  THE 
GMAC  MODEL  / 


I 
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No  missing 
values 


TEST  FOR  MISSING  VALUES 
AMONG  THE  COVARIATES  ^ 


Missing  Values 
among  co1  rariates 


REMOVE  MISSING  VALUES  FROM  THE  DESIGN  N 
MATRICES  OF  THE  VARIATE  MODEL,  REPLACE 
BY  ZEROES  AND  BUILD  EXTRA  COLUMNS  TO 
ACCOUNT  FOR  MISSING  VALUES  / 


j WRITE  THE  VARIATE  FORM  OF  THE 
MGMAC  MODEL  / 


COMPUTE  THE  VARIANCE-COVARIANCE 
MATRIX  FOR  THE  APPROPRIATE  VECTOF 
MODEL 


1. 


BUILD  THE  GROUPED  DATA 
\ FORM  OF  THE  MODEL 


/ 


WHITE  THE  GROUPED  DATA  j 
FORM  OF  THE  MODEL  / 


COMPUTE  THE  PARAMETER 
ESTIMATES 


z 


WRITE  THE  PARAMETER  ESTIMATES, 


READ  THE  NUMBER  OF  HYPOTHESES  TO  BE  TESTED 


READ  THE  TYPE  OF  HYPOTHESES  TO  BE  TAKEN 
EG.  FOR  ALL  VARIATES  OR  FOR  SIMPLE 
SPECIFIED  VARIATES 


READ  HYPOTHESIS  MATRIX 


/ 


/ WRITE  HYPOTHESIS  MATRIX,/' 


EXTEND  HYPOTHESIS  MATRIX  TO 
ACCOUNT  FOR  ADDITIONAL  PARA- 
METERS DUE  TO  MISSING  COVARIATES 


k 


CONSTRUCT  APPROPRIATE  HYPOTHESIS\ 
MATRIX  FOR  VECTOR  MODEL  OR  FOR 
SINGLE  SPECIFIED  VARIATE 


COMPUTE  THE  WALD  STATISTIC 


t 

v 


APPENDIX  B 


f\ 

i 

. 

PROGRAM  LISTING 

In  the  subsequent  pages,  a listing  of  the  FORTRAN  source  program 
is  given.  This  is  the  program  which  was  used  to  provide  the  printout 
for  the  sample  problem  in  Section  IV,  Volume  II. 

We  do  not  list  the  subroutines  used  from  the  Scientific  Subroutine 
Package  nor  the  International  Mathematical  and  Statistical  Library.  These  ‘ 

cure  available  in  the  appropriate  manuals. 


fortran  iv  ct 

0001 

0002 

0003 

000* 

0005 

0006 

000? 


0006 

0009 

0010 
001  1 
0012 
0013 
00!« 

0015 

0016 


0017 

0016 


0019 

0020 
0021 
0022 
0023 


RELEASE  2.0 

C 


RAIN 


OATE 


76320 


21/26/16 


C 

C 


REAL  *6  A( 100 .10 >.AJA( 10. 10) .AS( 100.10) .8(7.7) .BB( 100,  7 1 .BETA  ( 701 
*.8LM  100.2) .BSC  7 . 7 > , BT < 1 00 .6 1 .BSIGI7 .7) . BS IGB( 7 . 7 > ,BSI G1BI 7 . 7 1 . 
MBBS1GH7.7), BE  TAMXC26.7) 

INTEGER* A 1 2LK(  1001 • I TRT ( 1 0 0 1 . I R ( 7 1 . IRR ( 7 . 7 ) , IN ( 100 . 7 1 , I N0( 100 1 . 
XNCJ< 127 ) .NCL( 127) ,NCMK(9) .NCMC7.9) ,NCMC(  7,7,9) .NCTC( 7,7) ,NCT( 7) . 
8N*tA(  7)  , NMAC  (7,7  ) . MjM  R C « ( 9 1 
I N T E GER* 4 1 RON ( 2 6 > . I CCL( 26)  . IPRK00C7) ,K0RESP(7> 

RE AL * 4 FMT ( 20 ) .EPS,HYPID(5) 

REAL  *8  F(100.  10  >.V».D.SIGI.H(56.  70  > . HPRN  V(  56  . 70  > ,HPRH(  56.56)  . 
%H9ETA<56) .HBETAP(56> ,HP(8.  70 1 .PR002 ( 70 . 70 ) 

RE  AL  *8  R V ( 70 ) .RSS( 10. 10 ).RSSS( 100 ). S( 100. 10) ,T(  51B0) . SIG( 7. 7) . 
•SIGINV(7,7),TRT(100.4),U(100.100),V(100). 


REAL*8  XC100.3)  .YC100.7) ,YYS( 100 .7 ) . YS ( 1 00 .7 ) , Z(  1 00 ) . ZT< 1 00 ) . 
<ZY(100)«ZZ(7> 

c 

c ***  ARRAYS  RHICH  MAY  HAVE  TO  BE  RED  I MENS IONEO  TO  HANDLE  A LARGER  • 
c •**  OF  MISSING  VALUES  ARE  S A .A JA . AS , BET A . BET AMX .F.H .HPRNV , PR0D2. RV . 

C **•  RSS. RSSS.S, T ,NCKK .NCM, NCMC. I RCW. ICOL 

c 

C ***  CURRENT  DIMENSIONS  OF  ARRAYS)  INCREASE  THESE  ANO  THE  CORRESPONDING 

C ***  DIMENSIONS  TO  INCREASE  THE  CAPACITY  OF  THE  PROGRAM 

C 

NP  DIM  = 7 
NNO I M =*  1 0 0 
NBDIM  = 2 
NTDIM  = 4 
NKOI Ms  3 
NADIM  = 10 

NFOIM  = 10 
NPROD=70 
NHPRHD=S6 


EQUIVALENCE  ( T ( l ) , AS ( 1 ) ) , ( X( 1 ) ,S( 1 ) ) . CH( 1 ) ,U( 1 ) ) 

DATA  B/49*0.000/,ELK/200*0,000/.M/3920*0,000/. IN/700*!/, 
6PROD2/4  900*0  «0D  0/ ,RV/7O*O.CD0/.TRT/4OO*0, 000/ 

DATA  RSS/10040.0CO/ 


c 

NP 

EQUALS 

THE 

NUMBER 

OF 

OEPENOENT  VARIABLES 

c 

NT 

EQUALS 

THE 

NUMEER 

OF 

TREATMENTS 

c 

NO 

EQUALS 

THE 

NUMBER 

OF 

BLOCKS 

c 

c 

NK 

EQUALS 

THE 

NUMEER 

INDEPENDENT  VARIABLES 

C 

c 

READ 

INPUT  DATA 

READ (5,1 > NP .NT .NE.NK ,NN, IDGT .NMISS.EPS. 0. IPRKOD 
1 FORMAT ( 7 13. T 30. El  0.0.01 0.0, 71 2) 

VR 1 TE (6*990 1 ) 

9901  FORMAT! • 1* , ‘PARAMETER  VALUES  RE AO  FROM  FIRST  DATA  CARO  t •,////> 
WMITEI6, 9909 ) O 
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9909  FORMAT! *0* . • THE  CURRENT  VALUE  BEING  USED  FOR  THE  MISSING  CODE  IS  S 
<• .G16.6) 

WR I TE ! 6 , 99 10)  ICGT 

9910  FORM  AT ( • 0*  . • THE  VALUE  OF  IOGT  SUPPLIED  FOR  USE  IN  LPSOOR  SU0ROUTIN 
*E  IS  : * . IS  > 

WRITEI6.991  l ) F.FS 

9911  FORMAT! ‘O’ , ’THE  VALUE  OF  EPS  SUPPLIED  FOR  USE  IN  THE  DMFGR  SUBROUT 

FINE  IS  ! • , G 1 6 . 6 ) 

WHITEI6.9902 ) NMISS 

9902  FORMAT! *0* . *NUMEER  OF  MISSING  VALUES  IN  COVAR IATES  S'.IS) 
WR1TEI6.9903)  NFOIM.NP 

9903  FORMAT! 'O' . 'MAX  C IM  RESP  VECTOR  I • * 1 3 • • 5 CURRENT  DIM  RESP  VECTOR:*. 

• 13) 

WRITE ! 6.9904 ) NNOIM.NN 

990*  FORMAT! *0* . *MAX  NUMB  OBS : * . I 5 . • 5 CURR ENT  NUMB  OSSl'.IS) 
WRITEI6.9905)  NeOIM.NB 

9905  FORMAT!  *0*  . 'MAX  NUMB  BLOCKS ! * . I 5 , • 5 C URRENT  NUMB  BLOCKS:*. 15) 

WRITE  16 .9906 ) NTCIM.NT 

9906  FORMAT!  *0*  .' MAX  NUMB  T RT S i * . I 5 . • ! CUR RENT  NUMB  TRTSt'.IS) 
WRITF.t6.9907)  NKCIM.NK 

9907  FORMAT! • 0* .* MAX  NUMB  COVARS 5 •• IS ,*! CURRENT  NUMB  C0VARS:*.I5) 

WRITE  I b. 9900)  NACIM 

-9908  FORMAT  I *0 '.* MAX  NUMB  COLS  IN  MODIFIED  DESIGN  M ATR I X : * . I S ./ . /. I X . 

••CURRENT  NUMBER  CCLS  MAY  BE  SEEN  IN  LISTING  OF  MODIFIED  DESIGN  MAT 
URIX  TO  FOLLOW  LATER*,//) 

WRITEI6.9770 ) 

9770  FORMAT! *0' • 'PR! NT  OPTIONS  CHOSEN  FOR  THIS  PROGRAM:  0=N0PR I NT , 1 *PR t 

• NT*) 

9769  FORMAT! • 0* . 'PRINT  CPTICN  FOR  MAC  MODEL  : *.I5> 

9768  FORMAT! *0* . ‘PRINT  OPTION  FOR  GMAC  MODEL  i *.I5) 

9767  FORMAT! *0* . 'PRINT  OPTION  FOR  MGMAC  MOOEL  : *.I5> 

9766  FORMAT! • 0* . 'PRINT  CPTICN  FOR  SIGMA  t ITS  INVERSE  : *,I5) 

9765  FORMAT! *0* . 'PRINT  OPTION  FOR  MATRIX  MODIFIED  MOOEL  S * . 15 ) 

976*  FORMAT! *0* . ’PRINT  CPTICN  FOR  CEPENOENT  VARIABLES  AND  OESIGN  MATRIX 

• FOR  VARIOUS  MISSING  VALUE  PATTERNS  S »,I5) 

9763  FORMAT! *0*. ’PRINT  OPTION  FOR  BETA  VALUES  : *.I5) 

WRITEI6.9769)  I FRKOD ! 1 ) 

WR ITE! 6,9768 ) IFRK0012) 

WRITEI6.9767 ) IFRK00I3) 

• R I TE I 6. 9766 ) IFFKCO!*) 

WRITE!o.9765  ) I FRKOD ! 5 ) 

WRITFI6.976A ) IFF  KQD 16) 

WRITEI6.9763)  IFRK00I7) 

READIS.2)  FMT 

2 FORMAT ! 20 A* ) 

C WRITE  INPUT  DATA 

C . . 

WR  ITE!  6 • 3) 

3 FORMAT  1 1HI . 'LISTING  OF  INPUT  DATA’.///) 

IFINB.LE. 1 . AND.NT.LE.1 ) GO  TO  11 
IFINB.LE.l ) GO  TO  29 

IF1NT.LE.I ) GO  TC  28 

IKKK»0 

00  4 1*1 ,NN 

PEAOIS.FMT  ) IBLK!  I >.  ITRT!  I ) . 1 V!  I .J>  . J»1.NP ), ! X!  I ,J)  , J-t  ,NK) 

IKKK  = I KKK7 1 

WKITEIb.S)  IBLK! I ), ITRT! I ). < V! I , J) , J«I ,NP| ,( X! I . J) . J»I , NK) 


mms/m* 
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3 FORMAT ( 1 MO .214. 15F8. 2 ) 

C 

: SET  UP  BLOCK  DESIGN 

c 

00  7 11=1 ,NE 

IF( I BLK<  1 ) .EQ. I I > BLK ( I .1 1 >>1.000 
7 CONTINUE 

: SET  UP  TREATMENT  CESION 

00  6 1 1=1 , NT 

IF  { I TRT  ( n.EO.M)  TR  T ( I . I I > = 1.000 
6 CONTINUE 
A CONTINUE 

WR I TE ( 6. 9798  ) I KKK 

9798  FORMAT! • 0* .//•• C* . 13. • OBSERVATIONS  HAVE  BEEN  READ  FOR  THE  CURRENT 
* DATA  SET.  DOES  THIS  AGREE  WITH  THE  CURRENT  * OBS  GIVEN  EARLIER?* I 

BUILD  DESIGN  MATRIX  eY  AUGMENTING  BLOCKS.  TREATMENTS  AND  COVARIATES 

CALL  ARRAYI 2,NN,NB,NNDIM,NBDIM,BLK. BLK  > 

CALL  ARRAY  t2.NN.KT • NNO I M * NT  01M.TRT.TRTI 
CALL  CTIE(BLK.TRT.BT.KN.NB.O.O.NT) 

NBT  = N0fNT 
GO  TO  25 
11  I KKK  = 0 

00  41  1=1, NN 

RE  AL)  ( 5 , FMT  ) < Yt  I.J)  . J*1  .NP  ) . ( X ( I , J ) , J=  1 , NK  > 

IKKK= IKKK+1 

WRITE (6.95)  (Y<  I.J>,J=1.NP>.<X(I.J>. J»1,NK> 

95  FORMAT! 1H0. 15F8.2) 

41  CONTINUE 

WRITE(6.9798)  I KKK 
00  22  1=1, NN 
00  22  K= 1 , NK 
A!  I • K ) =X  ( I , K ) 

22  CONTINUE 
NM  = NK 
GO  TO  26 
29  I KKK  = 0 

DO  40  1=1, NN 

READ! 5. FMT ) I TR T ( I) , ( Y ( I , J > , J= 1 .NP). IX! I • J ) , J=1 .NK) 

1 KKK= IKKK4 1 

WRITE (6. 45 ) I TRT ( I ) . ( Y ( I . J ) . J> 1 . NP > . ( X ( t . J ) . J- 1 . NK ) 

43  FORMAT! IHO, 14, 13F8.2) 

DO  60  11=1, NT 

IF! I TRT! 1 ) .EO.I 1)  TRT! I, 1 1 )>1 .000 
60  CONTINUE 
40  CONTINUE 

WR I TE ( 6 , 9798 ) I KKK 

DO  24  11=1, NN 

00  24  K=  1 , NT 

QTI  1 I ,K>  = TRT(  II  ,K) 

24  CONTINUE 

CALL  ARRAY(2,KN,KT,NN0tM,NTDIM,BT,BT) 

NHT  » NT 
GO  TO  25 
28  !KKK=0 
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00  42  1=1 ,NN 

READ  ( 5.FMT  ) I BLK! I ) , ( Y< I , J)  , J = 1 ,NP) . <X< I . J >, J=l  ,NK) 
IKKK.=  IKKK+1 

WRITE (6 .46)  IBLK< I ) . ( Y< I. Jl . J=1 .NP) . (X! ! . J). J=1.NK ) 
46  FURMATI 1H0 , 1 4 . 15F8.2 ) 

DO  6 1 I I =1  . N0 

IFIIOLK!  I) .£0.1  I)  BLK(I.Il)  = l.ODO 
61  CJNT INUC 
42  CONTINUE 

WRITE16. 97981  I KKK 
00  27  1=1 .NN 
DO  27  K=l ,NB 
8T< I ,K)=BLK( l.K) 

27  CONTINUE 

CALL  ARRAY (2 .NN.NB.NNDXM. NbD XM.BT.8T > 

NOT  = NS 

25  CALL  A9RAY(2,NN.NK,NNDIM,NK0IM.X.X) 

CALL  CTIEIQT.X.A.M'i.NBT.O.O.NKl 
NM=NBT*NK 

CALL  ARRAY ( l . NN  . KM ,NNO I M . NA D 1 M . A . A 1 

N0M00=N8 
NT  MO  0=NT 

IF  (NB.EO.l)  NOI»CO=N0-t 
IF  (NT.EQ.l)  MTACC=NT-I 

NM  EQUALS  THE  NUMEER  OF  COLUMNS  IN  THE  DESIGN  MATRIX 
WRITE  THE  MATRIX  FCRM  OF  Y AND  A FOR  THE  MAC  MODEL 


26  WRITEI6.12) 

12  FORM  AT  11  HI ,• THE  VALUES  OF  Y ANO  A FOLLOW  MAC  MODEL*.///) 

IF  ( I PRX 00 (I  l.EC.O)  CO  TO  9797 

00  13  1=1. NN 

WRITEI6. 19)  I VI  I. J). J=1 ,NP) .( A( I ,J> . J=1 .NM) 

19  FORMAT! *0* . 1AF9. 2) 

13  CONTINUE 
GO  TO  9795 

9797  WR1TE<6.9796) 

9796  FORMAT! *0*.*  THE  ABOVE  LISTING  WAS  SUPPRESSED  1 IPRTXOO! 1 1*0* ) 
COMPUTE  THE  VARIATE  FORM  OF  V AND  A FOR  THE  GMAC  MODEL 
9795  00  6 J=1 ,NP 

WRITE  THE  VARIATe-WISE  FORM  OF  V AND  A FOR  GMAC  MODEL 
WR !TE!6. 161  J 

16  FORMAT  C 1 Ml . * THE  VALUES  OP  Y AND  A FOR  VAR  1 ATE  * . I 2 • • FOLLOW 
•GMAC  MODEL*.///) 

IF  ( IPRK0D(2).EC.C1  GO  TO  9794 
GO  TO  9792 
9T94  WR I TC < 6. 9793 1 

9793  FORM AT(*o'»*  The  ABOVE  LISTING  WAS  SUPPRESSED  i IPRTKOOI 2 )»0* > 
9792  !CT»0 

00  9 I>1  .NN 

IF(Y(  1.J1.E0.0)  GC  TO  9 

lCr»!CT*l 

Zf  1CT I a VI  I • J I 


63 


J&M  cof  ^ furnish®  to  Mfi  - 


FORTRAN  IV  fit  RELEASE  2.0 


DATE  » 76320 


21/26/16 


DO  10  K=1.NM 
F( ICT.K)  a A(I.K» 

10  CONTINUE 

IF  < IPRKODI2) .EO.C)  GO  TO  9 
WRITE(6.18)  Z(ICT),  (Ft  ICT.IC)  .X=l  . NM  ) 

18  FORMAT (• 0* . 14F9 .2 > 

9 CONTINUE 

NCT(J)  EQUAL  THE  NUMBER  OF  OBSERVATIONS  FOR  VARIATE  J , Jal  ..  . . .NP 


NC T ( J ) = ICT 
NBB=0 

IF(NMISS.EO.O)  CC  to  6000 
DO  3000  I - 1 . NN 

DO  3000  K = l.KWISS 

Brt< l .K ) = 0 . ODO 

3000  CONTINUE 


R fc"  Mu VE  MISSING  VALUES  FROM  A MATRICES  AND  REPLACE  BY  ZEROES 
BUILD  EXTRA  COLUMES  FOR  A MATRICES  TO  ACCOUNT  FOR  MISSING  VALUES 

KK=NHT*1 
DO  90  K^KK.NM 
L=NCT( J) 

NC  MI  J.K)=0 

DO  100  I a l.L 

IF(F( I.K1.NE.O)  GO  TO  100 

F ( I • K ) a 0. 000 

NCMIJ.K)  a NCM(J.K)  ♦ 1 
110  KXK=  NCM(J.K) 

B3( I .KKK)  3 1 .OCO 
NBB  = NBO  ♦ 1 
100  CONTINUE 
90  CONTINUE 

BUILD  NEW  A MATRICES  FOR  MGMAC  MODEL  BY  AUGMENTATION 
L a NCTI J 1 

nma(J)  eouals  the  number  cf  columns  in  the  design  matrix  for  variate  j 

NMA(J)  a NM  ♦ NEE 
do  1 30  KI>1 ,L 
NMP1  a NM  A 1 
NM A J a NMA(J) 

DO  130  K2  = NMPl.NMAJ 
NDUM  a K2  - NMP1  ♦ 1 
FIK1.K2)  a BQ(Kl.NOUM) 

130  CONTINUE 

WRITE  THE  VARIATE-WISE  FORM  OF  Y AND  A FOR  THE  MGMAC  MOOEL 


WR1TE(6.51)  J 

51  FORMATUHI  . 'THE  VALUES  OF  V AND  A FOR  VAR  IATE  • . 1 2.  • FOLLOW 
•MGMAC  MODEL*.///) 

IF  ( 1PRK0D(3).EC.C>  GO  TO  9791 
GO  TO  9790 
9791  WRITE (6. 9789) 
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0206 

9789 

FORM  A T ( • 0 • • • THE  ABOVE  LISTING 

WAS  SUPPRESSED 

0207 

GO  TO  6000 

02C8 

9790 

LL=NMA< J) 

0209 

DO  52  1=1 *L 

021  0 

MR  I T E ( 6 • 53 ) Z ( I ) • <F(  I.K).K=1.LL) 

021  1 

53 

FORM AT ( • 0» • 14F9.2 ) 

0212 

52 

CONT INUE 

0213 

6000 

NMA(J)  = NM  ♦ NEE 

0214 

L1=NCT( J ) 

0215 

Ml =NMA ( J ) 

21/26/16 


I PRTKOD ( 3 )=0 ' ) 


0216 
0217 
0 218 

0219 

0220 
0221 
0222 

0223 

0224 

0225 

0226 

0227 

0228 
0229 
0233 

0231 

0232 

0233 

0234 


C 

C COMPUTE  VARIANCE 

c 

C FINO  F'F 

C . 

C TEST  FOR  RANK  OF  f AND  COMPUTE  <F*F>  INVERSE=RINV 
C 

DO  3002  I = 1.M1 
DO  3002  L = liMl 
RS3<  I .L  ) = O.OOC 

DO  300 3 K = 1 .L  1 

IF(F<K, I l.EO.O.COO  OO.OR.F<K.L>.EO.O.OOO  00)  GO  TO  3003 
RSS(I.L)  = RSSM.L)  ♦ F(K.  I )*F(K,L> 

3003  CONTINUE 
3002  CONTINUE 

DO  900  I =1 , Ml 
DO  9 0 0 K = 1 . M 1 
IK1=< 1-1 >*M1  * K 

900  R5SS ( I K 1 ) =RSS< I .K ) 

CALL  DMFGR1 RSSS.K1.M1 .EPS, I RANK. I ROW, I COL) 

I R ( J ) = IRANK 
WRITE(6,S5)  J.IMJ) 

55  FORMAT ( I HO , • THE  RANK  OF  THE  OESIGN  MATRIX  FOR  VA RI ATE • , 1 5 , • FOR  TH 
*E  MG  MAC  MODEL  IJ'.IS.*  •«///> 

CALL  LPSDORIRSS  .V 1 . M 1 . NFD1M , RSS, IDGT.T, IER) 

WRITEI6.9864 ) IER 

9664  FORMAT! • 0* ./. • *,»ERRCR  CODE  FROM  LPSDOR  = *.16./) 


0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 


C 

C FORM  THE  PROOUCT  F(RINV) 

C 

9753  DO  3004  I = l.Ll 
DO  3004  L = 1 .Ml 
S(  I > L ) = 0.000 
00  3005  K = l.Ml 

IF(F ( I ,K ) .EO.O. CCO  00. CR.RSSIK. L ).E0. 0.000  00)  SO  TO  3005 
5(1. L)  = S(I.L)  ♦ F( I ,K)*RSS(K.L) 

3005  CONTINUE 
3004  CONTINUE 

c 

c 

C COMPUTE  FIRINVIF* 

C 


0243 

0244 

0245 

0246 
0 247 
0246 


00  3006  I s l.Ll 
00  3006  L = 1 .LI 
U(I.L)  s 0.000 
DO  3007  K = l.Ml 

IF(S( I.K). EO.O. COO  OO.CR.F(L.K) , EO.O. 000  00) 
U(I.L)  s U(I.L)  4 S( I ,K)*F(L.K) 


GO  TO  3007 


I 


) 


I 
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0253 

0254 


0255 

0256 

0257 

0258 

0259 

0260 
0261 


0262 

0263 

0264 

0265 

0266 


0267 

0268 

0269 

0270 

0271 

0272 

0273 
0 27A 

0275 

0276 

0277 

0278 
02  7 'y 
0280 
028  1 
0282 

0283 

0284 

0285 

0286 

0287 

0288 

0289 

0290 
0 29  1 

0292 

0293 

0294 
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3007  CONTINUE 
3006  CONTINUE 

C 

C COMPUTE  Z*IZ  = Z*Z 
C 

ZZCJ)  = 0.000 
00  320  I 1= 1 .L 1 
ZZ< J)  • ZZ( J)  ♦ Z( 1 I >**2 
320  CONTINUE 

c 

C COMPUTE  Z’FCRINV)F* 

c 

DO  3008  I .a  1 ,L1 
V( I ) = 0.000 

00  3009  L = 1 .LI 

1F(U(L.I I.EQ.O.COO  00)  60  TO  3009 
VC  I > = VC  I ) ♦ Z (L)»U(L,I ) 

3009  CONTINUE 

3008  CONTINUE 

c 

C COMPUTE  Z>FCRINV)F«Z 
C 

• 1 = 0.000 

00  3010  I a 1.L1 
« = W ♦ Z( I )*VC 1 > 

3010  CONTINUE 

SI6CJ.J)  = CZZCJ)  - WJ/CLl  - IR  C J ) ) 

c 

C COMPUTE  COVARIANCE 

c 

JJ  * J 1 

IFCJJ.6T.NP)  60  TO  S 
DO  41 1 J JJ  * JJ ,KP 
ICTC  = 0 
DO  4 09  I = 1 .NN 

IFCV ( I . JJJ) .EQ.C.CR.VC I.J1.E0.D)  60  TO  *09 
ICTC  » ICTC  * 1 
Z!  ICTC)  = V C I . J ) 

ZVCICTC)  = VII. JJJ) 

00  410  K a 1 . NM 
FC  ICTC.K)  a AC  I .K) 

*10  CONTINUE 
*09  CONTINUE 

NCTCCJ.JJJ)  » ICTC 
NBBC  = 0 

I FCNM ISS.EO.O ) CO  TO  6001 
00  3001  I a l.NN 
00  3001  K a l.NVISS 
BBC  I .K)  a 0.000 
3001  CONTINUE 

KK  a NOT  ♦ 1 
00  490  K = KK.KV 
L = NCTCCJ.JJJ) 

NCMCC J.JJJ.K)  a 0 
00  *100  I a I ,L 
IFCFC  I.M.NE.O)  GC  TO  4100 
F(  I.K)  > o.ODO 

NCMC ( J.JJJ.K)  a KC  VC C J.JJJ.K)  * 1 
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KKK  * NC  MC ( 

0296 

BOU.KKK)  = t.OCC 

0297 

NQQC  = NB3C  ♦ 1 

0 298 

4100 

CONTINUE 

0299 

490 

CONT INUE 

0300 

L * NCTC(J.JJJ) 

0 301 

NMAC ( J • JJJ i « N A ♦ N0BC 

0 30  2 

00  4 130  Ki  » l.L 

0303 

NMP1  = NM  4-  1 

0 30  4 

NMAJ  = NMAC(J.JJJ) 

0305 

00  4130  K2  = NMP1.NMAJ 

03C6 

NDUM  = K2  - N VP  I ♦ 1 

0307 

F(Ki«K2)  = BO(KItKDUM) 

0308 

4130 

CONT INUE 

0 309 

6001 

NMAC(J.JJJ)  a KV  ♦ N0BC 

0310 

L2  a NCTC(J.JJJ) 

031  1 

C 

M2  = NMAC(J.JJJ) 

C FIND  F'F 

C 

0312 

DO  3012  1 = l .M2 

0313 

00  3012  L = 1.M2 

0314 

RSSI I.L)  = O.OOC 

0315 

00  3013  K = 1 ,L2 

0316 

IF(F(K. I ).E0. 0. COO  00. OR 

0317 

RSS!  I.L)  = RSS!  l.L)  ♦ F!l 

0318 

3013 

CONTINUE 

0319 

3012 

C 

CONT INUE 

21/26/16 


TEST  FOR  RANK  OF  F ANO  COMPUTE  INVERSE 


DO  901  I = 1 .M2 
DO  901  K = 1 .M2 
IK  1 = ( 1-1  )*M2  ♦ K 
RSSS ( IK1 )=RSS! I .K) 

901  CONTINUE 

CALL  OMFGRIRSSS.M2.M2.EPS. (RANK. I ROW, I COL) 
IRR(J.JJJ)  > irank 

CALL  LPSOOR < RSS  , A2 . M2 . NFO IM , RSS . IOGT.T. IER) 
WRITE! 6.9864 ) IER 

FORM  THE  PRODUCT  F(RINV) 


DO  3014  1 » l ,L2 

DO  3014  L “ 1.M2 
S( I.L)  * 0.000 
DO  3019  K > 1.M2 

IF(F< I.Kl.EQ. 0.000  OO.OR.RSSCK.LI.EO.O.OOO  00)  CO  TO  3019 
S(I.L>  » S(I.L)  4 F ( I , K ) 4RSSIK.L > 

3015  CONTINUE 
3014  CONTINUE 

: COMPUTE  F(RINV)F» 


DO  3016  1 = 1 • L 2 
DO  3016  L * 1 ,L2 
U(  l.L)  « 0.000 
DO  J01 7 K « 1.M1 
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1F1S! I .KI.EO.O. 000  OO.OR.FIL.KI.EQ. 0.000  00)  GO  TO  3017 
U1I.L)  - U1I.L)  ♦ S( 1 .K)*F(L.K) 


3017  CONTINUE 
3016  CONTINUE 


COMPUTE  2 ’ I ZZ 


ZZ(J)  ^ 0.000 
DO  620  1 l>  1 ,L2 
ZZ(J)  » ZZ(J)  ♦ Z(II)*ZV(II) 
620  CONTINUE 

COMPUTE  Z*F!R!NV|F* 


00  3018  I a t , L 2 
V « I > X 0.000 
DO  3019  L = 1 ,L2 

IF(U<L. I l.EO.O.COO  00)  GO  TO  3019 
V < I ) = VII  ) ♦ Z ( L ) *U ( L . I ) 

3019  CONTINUE 
3018  CONTINUE 


C.  COMPUTE  Z*F|RINV)F*ZV 
W = 0.000 
OO  3020  I = 1 ,L2 
w * w ♦ ZV<  I ) • V < I ) 

3020  CONTINUE 

sigij.jjji  = (zz(j>  - *] 

SIG(JJJ.J)  = SIGIJ.JJJ) 
All  CONTINUE 
8 CONTINUE 


1RR1J.JJJ)) 


C *♦*  CALL  TO  SMOOTH  FOLLOWS  THIS  CARO 
C 

CALL  SMOOTH! S IG «NP« 1 ) 

C 

C WRITE  OUT  SIGMA  AAC  SIGMA  INVERSE 
C 

WRITEI6.69S) 

695  FORMAT (I  HI.  'THE  VALUE  OF  SIGMA  FOLLOWS1.///) 

IF  I IPRK0014 ) .£0.0)  GO  TO  9786 
DO  9788  JJJxI.NF 

WRITE(6. 9787 ) { SIG! JJJ.JJJJ) .JJJJ-l.NP) 

9787  FORMAT! • 0* . 10! 1 X.G12.6) ) 

9788  CONTINUE 
GO  TO  9784 

9786  WRITE16.9785) 

9785  FORM  AT ! * 0 1 * 1 ThE  ABOVE  LISTING  WAS  SUPPRESSED  1 IPRTKOOI4 >-0 • ) 
9784  CALL  LPSOORI S I G « AP. NP. NPO 1M. S 1GINV • 1 OGT . T* IER ) 

WR ITEI6.9864)  IER 
WRITE  16*699 1 

699  FORMAT ! 1 HI . ' THE  VALUE  OF  SIGMA  INVERSE  FOLLOWS1.///) 

IF  t IPRKOO ! 4 ) . EC. 0 ) GO  TO  9763 
OO  9782  JJJM.NF 

WRITE  I 6. 978 1 ) 1 SIGINVi JJJ.JJJJ) .JJJJ-l ,NP> 

9781  FORMAT! *0* . 101 1 X.G12.6) ) 

9782  CONTINUE 
GO  TO  9780 
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0 366 
0367 
0366 
0 369 
0 390 

0391 

0392 

0393 
039A 
0 395 
0 396 
0 397 
0396 

0399 

0400 

0401 

0402 

0403 


04C4 

0405 

0406 

0407 

0408 

0409 

0410 


04  11 

0412 

0413 

0414 

0415 

0416 

0417 
0416 
0419 


0420 

0421 

0422 


0423 

0424 


9763  MR I T E I 6 . 97  0 9 1 

C 

C REMOVE  MISSING  VALUES  FRCM  A ANC  REPLACE  BV  ZEROES 

C BUILD  EXTRA  COLUMAS  FOR  A MATRIX  TO  ACCOUNT  FOR  MISSING  VALUES 

C 

9780  NBB  s 0 

IFINMISS.EO.O)  GC  TO  6002 
00  5005  1 = 1 ,NK 

00  5005  K a 1 , AMISS 
B8( I »K)  « 0.000 
5005  CONTINUE 

KK  = NOT  4 1 

DO  1011  K a KK .AM 
NCMK1KI  * 0 
00  1012  I > l.NA 
IF( A ( I ,K I .NE .0 1 GO  TO  1012 
A 1 1 • K 1 a 0.003 
NCMK(K)  a NCMK1KI  4 1 
KKK  a NCMK1KI 
BB< I .KKK I a 1 . OCO 
NbD  a NOB  4 1 

1012  CONTINUE 
1011  CONTINUE 
C 

C 0U1L0  NEW  A MATRIX  FOR  MATRIX  MODIFIED  MOOEL 
C 

C NA  EQUALS  THE  NUMBER  OF  COLUMNS  IN  THE  MODIFIED  DESIGN  MATRIX 
C 

NA  = NM  ♦ N8Q 
00  1114  K 1 a 1 .AN 
NMP1  a NM  4 1 

DO  1114  K2  a NMFl.NA 
NDUM  a K2  - NMP1  4 1 

AIK1.K21  a BQIK1.A0UM) 

1114  CCNTINUE 
C 

C WRITE  THE  MATRIX  ACCIFIED  FORM  OF  V AND  A 
C 

NRITEI6. 1020 1 

1020  FORMAT! 1 HI  .' THE  VALUES  OF  V AND  A FOLLOW 
•MATRIX  MODIFIED  MODEL'.///! 

IF  I IPRKODI51.EC.01  GO  TO  9779 
00  1021  I a 1 . A A 

WRITE <6.  1022  1 (VII  . J 1 • J= 1 «NP I . I A! I .Kl.Kal.NAl 
1022  FORMAT! 'O' .14F9. 21 

1021  CONTINUE. 

WR ITCI6.97 29  I 

9729  FORMAT!  '1'.'  ••••••  DEPENDENT  VARIABLES  AND  DESIGN  MATRIX  FOR 

X THE  VARIOUS  GROUPS  CORRESPONDING  TO  DIFFERENT  PATTERNS  OF  MISSING 
S VALUES'  I 
NDIGHPa2**NP  - 1 
WR1TEI6. 97031  NCI  GRP . NDI GRP 

9703  FORMAT! 'O' . •••**  THERE  ARE  '.14.'  DIFFERENT  POSSIBLE  GROUPS  HUMMER 
XED  FROM  | TO ',14,'  THOUGH  IN  GCNERAL  NOT  ALL  GROUPS  WILL  APPEAR',/ 
M/.'O','***  WHIC4  ONES  OCCUR  OEPENOS  ON  THE  PATTERN  OF  MISSING  VALU 
<ES. HOWEVER  THE  TCTAL  » OBSERVATIONS  IN  ALL  CROUPS  MUST  EQUAL  NN • I 
GO  TO  6002 
9779  WRITEI6. 977SI 


J 
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0 425 

9778  FORMAT! *0*.*  TEE  ABOVE  LISTING  WAS 

C*?6 

WRITEI6.9729) 

0427 

6002  NA  s NM  ♦ NBO 

042ft 

DO  1019  J s 1 >NP 

0429 

B(J.J)  z 1.000 

043C 

1019  CONTINUE 

0431 

DO  1000  I = l.NN 

0432 

IND< I ) = 0 

0433 

00  1001  J = 1 * NP 

0434 

IF(Y< I . J) .NE.OI  GC  TO  1030 

0435 

INI  I . J)  z 0 

0436 

1030  L = NP  - J 

0437 

I NO  ( I ) = I NO I I ) ♦ IN(I.J)*2«*L 

0430 

1001  CONTINUE 

0439 

1000  CONTINUE 

0440 

XP  = 2*«NP  - 1 

044  l 

DO  1002  LL  * l.KP 

0A42 

NCL(LL)  = 0 

0443 

00  1003  1 « l.KA 

0444 

IF( I N0( I l.NE.LL ) GO  TO  1003 

0445 

NCL ( LL  > = NCLILLI  ♦ 1 

0446 

ICT  z NCLILLI 

0447 

DO  1004  J » l.NP 

0440 

VSI ICT.J)  = VI I .J) 

0449 

1004  CUNT  I NUE 

0450 

DO  1 005  K z 1 ,N A 

0451 

AS!  ICT. X)  = All  .X ) 

0452 

1005  CONTINUE 

0453 

1003  CONTINUE 

0454 

IFINCLILL) .EO.O  ) GO  TO  1002 

0455 

JCT  z 0 

0456 

DO  1 007  J z l.NP 

0457 

IFIVSIl.JI  .EO.O  GO  TO  1007 

0450 

JCT  z JCT  ♦ 1 

0459 

DO  1008  JJ  * l.AP 

0460 

bSIJJ.JCT)  z BIJJ.J) 

0462 

tooe  continue 

0462 

1007  CONTINUE 

0463 

NCJ(LL)  » JCT 

c 

C WRITE  THE  GROUPED  DATA  FORM  OF  V AND  A 

0464 

X z NCLILLI 

0465 

JJ  * NCJILLI 

0466 

WRITEI6. 1032)  LL.K.JJ 

0467 

1032  FORMAT! 1 HO. • THE  DEPENDENT  VARIABLES  < 
*X  FOLLOW  FOR  GRCLP * . IS . 1 WHICH  HAS*. 
* I ATES  * , /// I 

0468 

IF  I IPRK00I6) .EO.O)  GO  TO  9777 

0469 

DO  1 024  I z 1 ,X 

04  70 

WRITE! 6.  1025)  ! VSI I . J) .J>1 .NP) , 1 AS! I 

0471 

1023  FORMAT! *0* . 14F9  .2 ) 

0472 

1024  CONTINUE 

0473 

GO  TO  9776 

0474 

9777  WRITEI6.9773) 

0475 

9775  FORMAT! *0*.*  THE  ABOVE  LISTING  WAS 

C 

76320  21/26/16 


IPRTXODI  5 >*0  • ) 


FOR  THE  MODIFIED  MATRIX  MODEL 


I !PRTK00<6 )»0* I 


COMPUTE  THE  SUM  OVER  ALL  GROUPS  OF  I B< I NV I 6 • *S I GMA4B » I XA • I V 
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fortran  iv  at  uuui  1.0 


•I/2B/IB 


C COMPUTE  B'ASIGFA 

C 

9776  NAP  » NA4NP 
NJ  > NCJ(LL ) 

OO  1036  I * l.NJ 

00  1036  L * 1 • NF 

BSIG(I.L)  s o.OCO 
00  2066  K = 1 »NF 

IF(OS(Ki 1 ) .EQ.0.000  00)  GO  TO  2066 
BSIG(I.L)  » BSIG(I.L)  4 BS(  K . I ) *S IG< K.L > 

2066  CONTINUE 

1036  CONTINUE 

t COMPUTE  B**SIGVA*Q 

c 

DU  1037  I s l.NJ 

00  1 037  L * 1 .NJ 

BSIGDI 1 ,L)  - O.CDO 
00  2067  K » t . NF 

IFIBSIK.LI.EO. 0.000  00)  GO  TO  2067 
BSIG»(I,L)  - BSIGB(I.L)  4 BSIG( I .K)*B3(K.L) 

2067  CONTINUE 

1037  CONTINUE 

c 

: COMPUTE  t NV ( E 1 AS  I GMA4B) 

C 

IFINJ.EO.l)  GO  TO  1038 

CALL  LPSOORIBSIGB.NJ.NJ.NPOIM.BSIGB. IDGT.T.1ER) 
VRITEI6. 986 A ) IER 
GO  TO  1039 

103B  BSIGB(NJ.NJ)  » 1.000/BSIG8INJ , NJ ) 

C 

: COMPUTE  B ( 1 NV IB*  *S 1 GMA  4B ) ) 


1039  DO  1040  I a>  i.NF 
00  1040  L = 1 .NJ 
BBSIGKI  ,L>  * 0.000 
00  2070  K s 1 , N J 

IF1BS1 I .K> .EG. 0.000  00)  GO  TO  2070 
BBS1G111.L)  3 BBSIGI(I.L)  4 BS( 1 .K)*BSIGB(K.L) 
2070  CONTINUE 

1040  CONTINUE 


<B< INVIB'ASI GMA  46 ) ) XA  • ) 


NL  » NCL(LL) 


: GET  COLUMNS  OF  V 

00  1042  \ * l.NL 
00  1042  J » l.NJ 
TVS! 1 .J)  > O.ODC 
00  2072  K « 1 , NF 

IF1BS1K. J) .EQ. 0.000  00)  GO  TO  207S 
VTS(I.J)  3 VVSII.J)  4 VB(1.K)9RS(K.J) 
2072  CONTINUE 
1042  CONTINUE 
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0516 

0517 
05  8 
0 510 

0520 

0521 

0522 

0523 
052* 
0525 


0526 

0527 

0528 
0520 

0530 

0531 

0532 

0533 


053* 

0535 

0536 

0537 

0538 

0539 
05*0 
05*1 
05*2 


05*3 

05*4 

05*5 

05*6 

05*7 

05*8 

05*9 

0550 

0551 

0552 

0553 
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c 

C COMPUTE  B ( !NV (e **S IGMA*B> )XA» >V 

c 

11  = 0 

DO  1043  l = l.NP 
00  1 0*3  K = 1 ,N* 

11=11+1 
DO  10*3  J = 1,NJ 
00  2073  L = 1 . NL 

IF< BBSIGI ( I , J ) .EO.O.ODO  00 . OR . AS ( L • K ) . EO • 0 • ODO  00)  GO  TO  2073 
RV(II)  = R V ( II)  ♦ BBSIGI < I . J ) +ASIL.K >*VVS<L. J ( 

2073  CONTINUE 
10*3  CONTINUE 

c 

C COMPUTE  THE  SUM  CVER  ALL  GROUPS  OF  ( ( B< I NV( B* *S I GMA*B ) 8 • >XA • A) 

C AND  COMPUTE  THE  INVERSE 
C 

C COMPUTE  B( INV(B'*SIGMA*B)B' 

C 

DO  1 051  I » l.NP 
DO  1051  L = l.NP 
BSIG10I I .L ) = 0 .OCO 
DO  2081  K = l.NJ 

1FIBSIL.K) .EO.O.ODO  00. OR. BBSIGI I I .K ) .EO.O.ODO  00)  GO  TO  20BI 
USIGIO(l.L)  = BSIGiail.L)  + BBSIGI1 I.K)*BS(L.K) 

2081  CONTINUE 

1051  CONTINUE 
C 

C COMPUTE  A* A 

C 

NL  = NCL  ILL  ) 

DO  1052  1 = 1 . N A 
DO  1 052  L = l.NA 
A JA( I ,LI  = 0.000 
00  2082  K = l.NL 

IFt ASIA. I ) .EO.O .OCO  00)  GO  TO  2082 
AJA(I.L)  = AJA(l.L)  + AS(K.I)*ASIK.L) 

2082  CONTINUE 

1052  CONTINUE 
C 

C COMPUTE  I B( INV(8**SIGMA *B )6*)XAaA 

C 

C PUT  THE  ABOVE  INTC  A MATRIX 
C 

11  = 0 

DO  105*  I = l.NP 
DO  105*  K = l.NA 
KK  = 0 
11=11+1 
DO  1 05*  J = 1 . NP 
DO  105*  L * l.NA 
KK  = KK  + 1 

PH0D2 (II *KK ) = FPCD2 I I I.KK ) + BS IGIB(I.J) *AJA( L. K ) 

105*  CONTINUE 
1002  CONTINUE 
C 

C COMPUTE  THE  INVERSE  OF  THE  ABOVE  SIAM 
C 
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055 A 

0555 

0556 

0557 


0558 

0559 

0560 
056  1 

0562 

0563 
0 56A 

0565 

0566 

0567 

0568 

0569 


0570 

0571 

0572 

0573 


0S7A 

0575 

0576 

0577 


0578 

0579 

0580 

0581 

0582 

0583 

058A 

0585 

0586 

0587 


IV  G1  RELEASE  2.0  VAIN  DATE  ■ 76320  21/26/16 

CALL  LPSDOR!PROC2,NAP,NAP.NPRCO,PROD2, IDGT.T ,IE«> 

WRITE (6.9064)  IEP 
DO  1056  1 = l.NAP 
BETA ( I)  = 0.0  DO 
C 

C COMPUTE  BETA  ESTIVATES 
C 

DO  1056  L = l.NAP 

BETA ( 1 ) = BET  A ( I ) ♦ PR0D2I I ,L )*RV!L) 

1056  CONTINUE 

DO  9720  18=1, NP 
DO  9728  19=1 .NA 
I 7= ( 10-1  ) ANA* I 9 
9728  BETAMX! 19. 10 >=BETA( 1 7) 

WRITEI6.1067) 

1067  FORMAT  ! • l 1 , • ******  LISTING  OF  PARAMETER  ESTIMATES  IN  THE  ORO 

#£R  : BLOCKS,  TR  T S . REGRESSION  COEFFICIENTS  OOWN  THE  PAGE  *•••**•) 
IF  < lf>RKOO(7).EC.O)  GC  TO  9774 
WR  I T E ( 6 , 1062) 

1062  FORMATI'O*.*  VAR*1  VAR»2  VAR#  3 

< VAR*4  VARK5  VAR«6 

X VAR #7  * ) 

00  1 060  10=1  ,NM 

WRITE (6.9727)  { eETAMX! 18. 19) . 19=1 , NP ) 

9727  FORMAT (,0,»7( 1X.G18.6)) 

1060  CONTINUE 
C WRITEI6. 1926) 

C 1926  FORMAT!//, • ******  ESTIMATES  OF  DUMMY  PARAMETERS  1NT 

C XRODUCED  BECAUSE  OF  MISSING  VALUES  ****«*•) 

C IF  (NM.EQ.NA)  GC  TO  19)9 

C NMP1=NM*1 

C DO  1925  I8=KMP1 .NA 

C WRITE(6.9727 ) ( EETAMXI 18.19) « 19= 1 , NP ) 

C 1925  CONTINUE 
C GO  TO  1917 

C 19 1 9 WR1TE16, 1910) 

C 1918  FORMAT  I * 0 # * • THERE  ARE  NO  DUMMY  PARAMETERS  ! THERE  WER 

C #£  NO  MISSING  VALUES  )•) 

GO  TO  1917 
977*  WRITEI6.9773) 

9773  FORMAT! •O'.*  THE  ABOVE  LISTING  WAS  SUPPRESSED  S IPRTKOOI 7 )*Q • | 
1917  CONTINUE 
C 

C READ  IN  THE  SPECIFIED  HYPOTHESES  MATRICES 
C 

WR I TE ( 6 » 9756 ) 

9756  FORMAT! ******  HYPOTHESIS  TESTING  SEC 
ST I ON  ******  •,///) 

READ ( 5 , 2020 ) AUVHVP 
2020  FORMAT! 15) 

WR I TE ( 6 , 9772 ) NUVhYP 

9772  FORMAT! *0* . 15. • HVPOTHESIS  MATRICES  SHALL  BE  USED, IN  TURN. FOR  COMP 
•UTING  CHI  SQUARE  STATISTICS*.///) 

WR  I TE ! 6 , 9771  ) A* 

9771  FORMAT!  *0*  . 'F.AC*-  MATRIX  SHOULD  HAVE  *.15,*  COLUMNS  AS  FOLLOWS  ) •) 
WR I TE ( 6 .9759 ) NEMCO 

9759  FORMAT! *0*.*  THE  FIRST  * . IS . • COLUMNS  CORRESPOND  TO  BLOCK  PARAMITE 
MRS  * ) 
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0588 

0589 

0590 

0591 

0592 

0593 
0 594 

0595 

0596 

0597 

0598 

0599 

0600 
06  01 
0602 

0603 

0604 


0605 

0606 
0607 
0 60S 
0609 
06  10 


061 1 
0612 
0613 
06  14 

0615 

0616 
0617 
061  8 

0619 

0620 
0621 
0622 


! 


0623 

0624 

0625 


VRITEI6* 97 5 E ) NTMOD 

9758  FORMAT! • •••  THE  NEXT  *.15.*  COLUMNS  CORRESPOND  TO  TREATMENT  PARA 
XMETERS*  > 

WRITEI6, 97571  NK 

9757  FORMAT!*  *.*  THE  LAST  *,I5«*  COLUMNS  CORRESPOND  TO  COVARIATE  COEF 
XFICIENTS* ) 

DO  2000  I = 1 .NLMHYP 

READ ! 5 .202 1 1 NUMRCMl I > .HYP I O. KOOOVL • ! K OR ESP! 191. I9»l ,NP> 

2021  FORMAT! 15. 5A4.2CI21 
NR  * NUMR0AII1 
VR1TEI6, 1999)  I.HVPIO 

1999  FORMAT! 1H1 . *****  LISTING  OF  HYPOTHESIS  MATRIX*. 15.*  FOLLOWEO  BY  IT 
»S  EXTENSION  FOR  MODIFIED  MODEL . HYPO T H ID  IS!  *.5A4,///) 

DO  2001  J = l.NR 

RLAOI5.2C 22)  (H ( J ,K ) ,K=t ,NM ) 

2022  FOMMAT120F4 .2 ) 

FRITE16. 1998)  IH| J.K) ,K*1 ,NM) 
l 99!)  TORMATl  ■ O*  . 20  I 1 X.F5. 2 ) ) 

2001  CUNTINUF. 

NHP  * NR*NP 

C 

C NRP=NR*NF  ANO  NAp*N A*KP 

C- 

C EXTENO  HYPOTHESIS  MATRIX  TO  ACCOUNT  FOR  AOOITIONAL  PARAMETERS 
C RESULTING  FROM  MISSING  INDEPENDENT  VARIABLES 
C 

IFINM.E0.NA)  GO  TC  2002 
DO  2003  J * l.NR 

00  2003  K * 1 .Nee 

KK  : NM  » K 

HIJ.KKI  * O.ODO 

2003  CONTINUE 

C 

C BUILD  HYPOTHESIS  MATRIX  FOR  PARAMETER  VECTOR 
C 

2002  DO  2004  L = I.NF 

00  2004  J * l.NR 

JJ  * J ♦ NR*!L  - 1) 

00  2004  K * 1 • N A 
KK  * K ♦ NA*!L  - 1) 

HIJJ.KK)  * HI J.K) 

2004  CONT INUE 

DO  4000  J*I ,NRP 

MR ITEI 6.9755 ) J . ! H I J . K ) , K»1 . N AP ) 

9755  FORM  ATI • *,//.*  * . • ROW! * . 1 2 . • ) 1 * . 20 ( 1 X . FS. 2 ) ) 

4000  CONTINUE 

IF  (KOOOVL. EQ.O)  GO  TO  9749 

C 

C COMPUTE  THE  TEST  STATISTIC 

C 

C COMPUTE  H*  APR  I NV 

c 

00  20 05  J - I.NRP 
00  2005  L * 1 .NAP 
HPRnvIJ.LI  ■ O.CCO 
DO  2006  K * 1 .NAP 

IFfHIJ.Kt. EG. O.CCO  OO.OR.PR002!K.L).EO.O.ODO  001  60  TO  S006 
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0629 

0630 


0631 
06J2 
0633 
063* 
0*  35 

0636 

0637 
0636 


0639 

06*0 


06*  1 
06*2 
06*3 
06** 
0 6*5 
06*6 
0647 


06*6 

06*9 

0650 

065 1 

0652 


0653 

065* 

0655 

0656 

0657 

0658 

0659 

0660 


0661 

0662 

0663 

066* 


2006  CONTINUE 
2005  CONTINUE 

c 

C COMPUTE  M'*PRINV«h 

c 

00  2007  J » 1 « NAP 
OO  2007  L » l.NRP 
HPRHU.Ll  * o.occ 
00  2008  K s 1 . NAP 

I F I NPRNtf (J.K).EQ.C.OOO  OO.OP.HIL.Kl.EQ.O.OOO  00)  CO  TO  2008 
MPRHIJ.L)  ■»  HPRMJ.LI  * HPRNV < J.K) *HIL .K ) 

2008  CONTINUE 

2007  CONTINUE 
C 

C COMPUTE  THE  INVERSE  OF  H>*PRINV*H 

c 

CAU.  LPSDORlHPRt-.KRP.NRP.NHPRHO.HPRH,  IDGT..T,  1£R  ) 

MRITE(6.9864 > IER 
C 

C COMPUTE  M**OETA 
C 

00  2009  J = 1 .NAP 
MUETA(J)  = 0.000 
00  2010  L * l.NAP 

IF(H(  J.D.EO.O.CDO  00)  GO  TO  2010 
HBETA(J)  = HBETAIJ)  ♦ HI J ,L ) *EETA(L) 

2010  CONTINUE 

2009  CONTINUE 

C 

C COMPUTE  CH**BETA) •( INV I H • *PR I NV *H ) ) 

C 

00  201 1 J » l.NRP 
HUETAPIJ)  * O.OCO 
00  2011  L = l.NRP 

HBETAPIJ)  » HOE  TAPI J ) ♦ H0E  TAIL) *HPR  H I J , L ) 

2011  CONTINUE 
C 

C COMPUTE  (H»BET A)* I 1NV« h»*PRINV*H) ) |H*8ETA) 

C 

■ALO  « 0.000 
00  2012  J 3 i ,NFP 

• ALO  > WALO  * H EE  TAPI J ) 4HBET  A t J! 

2012  CONTINUE 

MR  1 TE 1 6 • 20 13)  I .NP .WALO 

2013  FORMAT  I 1 HO *//.lNO»*THE  WALD  STATISTIC  FOR  HYPOTHEC  IS • . I* . • FOR  ALL 

• •»!*,•  RESPONSE  VARIABLES  SIMULTANEOUSLY  IS  i '.FIS.*.//) 

WRITE  16*2015 ) NR  F 

2018  F0RMAT1 *0* . • ITS  ASYMPTOTIC  DISTRIBUTION  UNDER  THE  NULL  HYPOTHESIS 

• IS  CHI-SQUARE  WITH  *.IA,»  DEGREES  OF  FREEOOM* .///) 

C 

C COMPUTE  THE  TEST  STATISTIC  FOR  THE  TESTS  ON  INDIVIDUAL  RESPONSE  VARIABLES 

C 

97*9  WRITEI6.9750)  I.HYP1D 

9750  FORMATIM'.*  •♦•*••  RESULTS  OF  HYPOTHESIS  TESTS  ON  INDIVIDUAL  VARl 
BATES  FOR  HYPOTM  FATRIX'.IS.'  WITH  IDS  * • 5 AA . * ****** • . /// I 
WRITE (6. 9708 ) NF. IKORESPI 19) . I9«l . NP ) 

9708  FORMATI//.'  •••  THE  OPTIONS  FOR  THE  ».I3,*  INDIVIDUAL  VARIATE 

*9  ARE!  • .1013./) 
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0665 

0666 


0667 

0668 

0669 

0670 


0671 

0672 

0673 
067* 

0675 

0676 

0677 

0678 


0679 

0680 
0681 
0682 

0683 

0684 
0635 
0686 


0687 

0668 


0689 

0690 

0691 

0692 

0693 

0694 
0693 


0696 

0697 

0698 

0699 

0700 


0701 


DO  9748  IV®1 .NR 

IF  (KDRESPI IVJ.EQ.O)  GO  TO  9748 
C 

C EXTRACT  APPROPRIATE  PART  OF  THE  H MATRIX  CONSTRUCTED  EARLIER 

C 

DO  9747  IVl =1 .NR 
00  9747  IV2*I.NAP 
1V3=! IV-1 l*NPtl VI 
9747  HP ( I VI . I V2 1— HI  I V3.1V21 

C 

C COMPUTE  H • *PR I N V 

c 

00  9 746  J®1.NR 
00  9746  L® 1 .NAP 
HPRNVI J.L1 = 0.D0 
DO  9745  K® 1 ,NAP 

IF  (HP(J.K) .EO.O. CO. CR.PP002IK.L l.EO. 0.001  GO  TO  9745 
HPRNVI J.L)=HPRNVC J.LI  ♦ HP( J . K ) 4PR0D2 ( K . L > 

9745  CONTINUE 

9746  CONTINUE 

c 

C COMPUTE  H**PRINV*F 

c 

00  9744  J® 1 . NR 
00  974*  L=l .NP 
HPRHI J .L ) *0 • 00 
00  9743  K® 1 . NAP 

IF  (HPRNVI J ,K). EC. 0.00. OR. HP(L.K).EO. 0.00)  GO  TO  9743 
HPRHI  J.L>=HPRH(  J.L)  4 HPRNVI J.K) 4HPIL.K) 

97*3  CONTINUE 
9744  CONTINUE 

C 

C COMPUTE  THE  INVERSE  OF  H'APRINVAH 
C 

CALL  LPSDORIHPRF.NR.NR.NHPRHO.HPRH. IDGT.T.IER) 

WRITE (6 .9864 ) IER 

c 

C COMPUTE  H'*3ETA 
C 

DO  9742  J®1.NR 
H3ETA( J>*0.00 
DO  9741  L® 1 .NAP 

IF  ( HP( J.L) .EO. 0.00)  GO  TO  97*1 
HbETAI  Jl-HBETM  J)  4 HP  ( J , L ) *BET  A <L  ) 

97*1  CONTINUE 
9742  CONTINUE 
C 

C COMPUTE  (H**8ETA)  •(  !KV(M**PRINV*H) ) 

c 

00  9740  J»1 »NR 
HbET AP( J)»0.00 
00  9740  L«l .NR 

H0ETAP1 J)=HBETAF( J)  4 H8ETA(L)*HPRH( J.L) 

9740  CONTINUE 

C 

C COMPUTE  (H*BETA ) • ( I N V • H • *PR IN V*H ) MH*8ETA) 

c 

■ALOsO.OO 
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070* 
0703 

0706 

0707 

0708 

0709 
071  0 

0711 

0712 

0713 
071* 
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DO  9739  J»I .NR 

9739  VALO=»ALO  ♦ MBS TAR { J ) *HBETA ( J 1 
WRITEI6. 97381  t, IV, HALO 

9738  FORMAT! *0* .* THE  »ALO  STATISTIC  FOR  HYPOTHESIS  •.!*,•  RESTRICTED  TO 
X RESPONSE  VARIATE  •,!*.•  ONLY  IS  ! '.FIS.*,/) 

MR I TE ( 6 , 9737 1 NR 

9737  FORMAT! • 0* ,• ITS  ASYMPTOTIC  DISTRIBUTION  UNDER  THE  NULL  HYPOTHESIS 
* IS  CHI-SQUARE  »ITH  •,!*,•  DEGREES  OF  FREEDOM* ,////) 

97*8  CONTINUE 

DO  9736  19*1, NHRRHD 
DO  9736  18=1,NPR00 
9736  HI  19 . 18  1=0.00 
2000  CONTINUE 
9999  STOP 

DEBUG  SUBCHK 
END 


FOB 
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MAIN 


DATE  * 76320 


21/26/16 


•OPTIONS  IN  EFFECT*  NOTERM, ID.E0COI C . SOURCE. NOL 1ST. NODECK .LOAO. NOMAP. NOTEST 
•OPTIONS  IN  EFFECT*  NAME  » MAIN  , LINECNT  = 60 

•STATISTICS*  SOURCE  STATEMENTS  * 715. PROGRAM  SIZE  * 376500 

•STATISTICS*  NO  DIAGNOSTICS  GENERATED 


78 


~ - 





MfsTc  MGR  IS  BEST  QUA LIT!  JBACTlWBLB 


FORTRAN  IV  Cl 


RELEASE  2.0  SMOOTH 


DATE  - 7*320  21226/1* 


0001 

C 

c •*• 
c **• 
c *** 
c *** 
c *** 
c 

c *** 
c *** 
c **• 
c *«• 
c *** 
c *** 
c *** 
c *** 
c *** 
c *** 
c 


SUBROUTINE  SMOOThlAI.N.KOOl 

THIS  SUBROUTINE  ACCEPTS  A DOUBLE  PRECISION  INPUT  MATRIX  AI  OF  DIMENSION 
N BY  N AND  CHECKS  TO  SEE  IF  IT  IS  POSITIVE  SEM IOEF IN  I TE . IF  SO  IT  IS  LEFT 

intact  and  retu-neo.  if  not  it  is  smoothed  by  the  least  squares  algorithm 

I NOT  GENERALIZED  LEAST  SQUARES  ) OF  SCHWERTMAN  ANO  ALLEN  (UNIV  KENTUCKY 
DEPT.  STATISTICS  TECH  REPORT  *56 .SEPTEMBER  1973)  TO  THE  NEAREST  POSITIVE 
SEMIDEFINITE  MATRIX  AND  RETURNED. 

AN  OPTION  IS  AVAILABLE  TO  PRINT  BOTH  THE  MATRIX  SUBMITTED  TO  THE  ROUTINE 
AND  THE  ONE  RETURNED  BY  IT  I 

KCC  = 0 BYPASSES  THE  PRINT  OPTION 
KCD  = 1 INVOKES  THE  PRINT  OPTION. 


THE  PROGRAM  CAN  EE  MADE  TO  HANDLE  LARGER  MATRICES  BY  INCREASING  NDIM.  THE 
DIMENSIONS  OF  VA.  EG.  AND  AI  IN  THE  SUBROUTINE  SMOOTH.  THESE  SHOULD  MATCH 
THE  O l MENS  I ON  OF  THE  AREA  HOLDING  THE  INPUT  MATRIX  IN  THE  CALLING  PROGRAM. 


0002 

REALMS  VA<7,7),TRC.EGC7.7),AU7.7) 

0003 

ND I M = 7 

0004 

DO  14  1=1 .N 

OOOS 

. 

DU  14  J=1,N 

0006 

14 

VA ( I ,J)=A1 ( I . J) 

0007 

IF  (KOO.EO.OI  GC  TO  13 

OOOS 

CALL  MATOUT  ( VA  .N.N. • INPUT  MATRIX  TO  SMOOTH  •.NDIM.NOIM) 

0009 

13 

TRC=0 . 0D0 

0010 

CALL  HOI AG(  VA.N.C.EG.NR.NDIM) 

001  1 

DO  15  1=1, N 

0012 

IF  < VA( l . I > .LT.C.OOOI  TRC»TRC4VA<t.n 

0013 

IS 

CONTINUE 

0014 

DO  16  1=1. N 

0015 

IF  (VA(I.I).LT.C.OOO)  GO  TO  17 

0016  ' 

16 

CONTINUE 

0017 

GO  TO  40 

00  1 B 

17 

DO  21  1=1, N 

0019 

DO  21  J= 1 , N 

0020 

21 

AI ( I , J)=0.000 

0021 

DO  20  1=1, N 

0022 

IF  < VA< I , I > .LT.0.000 > GO  TO  20 

0023 

00  19  J*1,N 

0024 

00  19  L=1.N 

0025 

AI  ( J ,L)=AI ( J.L)  4 EG(J.t)  * VA(t.I)  * EG(L.I) 

002* 

19 

CONTINUE 

0027 

20 

CONT I NUE 

0020 

40 

IF  (KOD.EO.O)  GC  TO  41 

0020 

CALL  MATOUTIAI. N.N. 'MATRIX  OUTPUT  BY  SMOOTH  '.NOIM.NOIMI 

0030 

41 

RETURN 

00  31 

OCBUG  SUBCHK 

0032 

END 
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0001 

SUBROUTINE  HOI  AG  ( H. N . IEGEN . U.NR , NN ) 

0002 

implicit  real*8  (a-h.c-z) 

0003 

01  MENS  ION  H(NN.NN)  , U(NN.NN)  .XII 00)  .10(100) 

0004 

IF(  I EGEN )50. 10.50 

0005 

10 

0040 1=1 .N 

0006 

0040 J=1 .N 

0007 

IF( I-J130.20.30 

0008 

20 

U(  I . J > = 1 .000 

0009 

G0TQ40 

0010 

30 

U< I • J ) =0.000 

00  1 1 

40 

CONTINUE 

0012 

50 

NR=0 

0013 

IF ( N— 1 ) 4 70 .470,60 

00  14 

60 

NMI 1 =N- 1 

0015 

0090  1=1 , NMI  1 

00  16 

X( 1 >=0.000 

0017 

IPLl=I+l 

0018 

0080 J=IPL1 , N 

0019 

IF(X< 1 )-0AaS(H( I ,J) >) 70.70.80 

0020 

70 

X( I ) =DABS( H( I.J1) 

0021 

I Q ( I ) = J 

0022 

80 

CONT INUE 

0023 

. 

RAP=7.450580596C-9 

0024 

HOTE  ST  = 1 .0038 

0025 

90 

001201=1 .NMI 1 

0026 

IF( 1-1 > I 10.  1 10, 100 

0027 

100 

IF ( XMAX—  X(  I )) 1 1 C. 120.120 

0028 

no 

XMAX=X< I ) 

0029 

IPIV=I 

0030 

JP I V = I 0(  I ) 

0031 

120 

CONTINUE 

0032 

IF (XMAX >470.470. 130 

0033 

130 

IF(HOTEST)  150, ISC. 140 

0034 

140 

IF(XMAX-HDTEST) 150.150.180 

0035 

150 

HO  1 M I N=D  ABS  ( H ( 1 .1  )) 

0036 

001701=2, N 

0037 

IF(HOIMIn-0ABS  (H( I . I > ) ) 170. 170, 160 

00  38 

160 

HO  I M I N=D  At!  S (H(  !,!>). 

0039 

170 

CONT  INUE 

0040 

HOTE  S T =HD IMIN4RAF 

004  1 

lF(HOTEST-XMAX)  18C, 470,470 

0042 

180 

NR=NR+1 

0043 

1F(H(  IP  I V.  IPI V)-M JPIV.JPIV ) >200 .190.200 

0044 

190 

S=1.000 

0045 

TANG  = OSIGN  (2.0C0.S) *H( IP I V . JP I V )/{ DABS ( H( IP  IV. IPIV )-H( JPI V, JPI V) 1 
140SQRT  ((H(IP1V.IF1V)-H( JPI V, JPI V) )**2+4.0DO»H( I P I V . JP I V ) ** 2 > ) 

0046 

GO  TO  210 

0047 

200 

T A NG  = D5 1 GN  ( 2.0C0. (H(  I P 1 V . IP  I V )-H( JP I V , JP I V ) ) )»H( I P I V, JP I V ) /< DABS! 
1H(  IPI V, IPI V)-H(  JPIV.JPIV)  )40SORT  ((H( IPIV, IPIV)-H( JPIV.JPIV ) )«*244 
2.0D0*H(  IPIV, JPI V >442)  ) 

0048 

210 

C0SINE=1 .O/DSQRT  (1 .0004TANG442 ) 

0049 

SINE=TANG*COSINE 

0050 

HI I=H( IPIV, IPIV ) 

0051 

H(  IPIV. IPI V)  = C05INE**2*(HII4TANG*(2.000*H< IP I V, JPI V ) 4TANG*H( JP I V, J 
IPIV)  ) ) 

0052 

hijpiv.jpiv)=cosine*42*( h( jpi v. jpi v)-tang*( 2. odo*h{ ipi v, jpi v> -tang 

1 * H I I ) ) 

0053 

h(ipiv,jpiv>=o.coc 

81 
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0054 

1F(H( IP  IV. IPI V) -HI JPIV, JPIV) ) 220.230.230 

OC  55 

220 

HTEMP=H(  IPI V.  IP  IV ) 

0056 

H<  IP  I V, IPI V)=H(  JF I V. JPIV) 

0057 

H(  JP  IV.JPIV  ) =HT  EMP 

0056 

HTEMP=DSIGN  (1.CCO.-SINEI4COSINE 

0059 

C0SINE=DABS  (SINE) 

0060 

SI NE=HTEMP 

0061 

230 

CONT l NUE 

0062 

003101=1 .NMI1 

0063 

IF! I- IPI V) 250.3 10 .240 

0064 

240 

IF! I-JP IV) 250. 2 10. 250 

006  b 

250 

1F(IO(I>-IPIV>2CC.270.260 

0066 

260 

IF! IQ ( I )-JP IV 13 10. 27 0.3 10 

0067 

270 

K=  IQ  I 1 ) 

0066 

280 

MTEMP  = I((  1 ,K  > 

0069 

HI  I , K ) =0.000 

0070 

I PL  1 * I ♦! 

0071 

XI  1 I =0  .000 

0072 

D0300J=IPLI «N 

0073 

IF! XI I l-OABS  IHII.JII 1290.290.300 

0074 

290 

XII) =DABS  (HII.J)) 

0075 

10(1  >=J 

0076 

■300 

CONT INUE 

0077 

HI  I . K ) =HTEMP 

0078 

310 

CONT I NUE 

0079 

xi  ip  i v > = o.ooo 

0080 

XI JPIV)=0.000 

0081 

00440  1=1  .N 

0062 

IF! I — IPIVI320.440.360 

0083 

320 

HTEMP=H( I . IPI V) 

0084 

HI  I , IP  I V )=COS INE4HTEMP  + SINE4HI I . JPI V ) 

0085 

I F I X I I l-OABS  (HII.IPIvni330.340.340 

0086 

330 

XI I ) =DABS  (HI  1 , IPIV1) 

0087 

1 Q I I )= IPIV 

0088 

340 

HI  I , JPI V >=-S INE*HTEMP  + CCSINE4M( I . JP 1 V 1 

0089 

IF  I XI  I ) -DA  IS  (H(  I , JP I VI)  >350.44  0.44  0 

0090 

350 

XI  1 ) = DAOS  I HI  I , JPIV I ) 

0091 

IQI 1 >=JP!V 

0092 

G0TQ440 

009  3 

360 

IF! I-JP! VI 370. 440. 400 

0094 

370 

HTEMP=M( IPI V. 1 ) 

0095 

HI  IP  IV. I )=COSINE*hTEMP*SIN£*H( I,  JPIV) 

0096 

IFIXI  IPIV)— OABS  (HI  IP  IV. I ) ) 1380.390.390 

0097 

380 

XI IPI V)=OABS(H( IPIV. I )) 

0098 

IQI IPIV>=! 

0099 

390 

HI  I . JPIV ) =-SINE»HTEMP*CCSINE*M( 1 .JPIV) 

0100 

IFIXI  I) -DABS! HI  I. JPIV) )> 390. 440. 440 

0101 

4 00 

HTEMP=H( IPIV. t ) 

0102 

hup iv. i )=cosine*htenp4sine*hi  jpiv.i  > 

0103 

IFIXI IPIV) -OABS  (HI IPIV, 1 )) 1410.420.420 

0104 

410 

X(IPIV)=DAQS  (HIIFIV.t)) 

0105 

IQI 1PIV)=I 

0106 

420 

HI JP I V. 1 )=-SINE*HTEMP>C0SINE4H( JPIV, I ) 

0107 

(FIX! JPI VI -OABS  (HI JPIV. I II >430,440,440 

0108 

4 30 

XI JP I V 1 =OAB  S (H(JPIV.I)I 

0109 

IQI JPIV)*! 

0110 

440 

CONT INUE 

0 1 1 1 

IF ( I EGEN )90 .450 ,90 

21/26/16 
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0112 

450 

00460 1=1 ,N 

0113 

HTEMP=UC I » 1 P I V ) 

Oil* 

U<  I.  1 PI  V ) =COS  1 NE*HTE  NP+S I NE*U ( I . JPIV1 

0116 

460 

U< I.JPIV)=-S I NE*HTEMP. COSINE *U1 I . JPI V ) 

0116 

GOT090 

0117 

470 

RETURN 

01IS 

DEBUG  SUBCHK 

01  19 

end 

21/26/16 
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0001 

0002 

0003 

0004 

0005 

0006 
0007 

oooa 

0009 

0010 
0011 


SUBROUTINE  MATOLT (A.M.N. TITLE. NR .MC ) 
RE  AL  * 8 A 

DIMENSION  AINR.UCI.TI TLE  C 6) 

1 FORMAT (,0,/,0* . 6 A4  ) 

2 FORMAT! • • «7G1 8.8) 

WRITE(6.1)  TITLE 
00  3 1=1 ,M 

3 WRITE(6.2>  ( A ( I , J ) , J= 1 ,N ) 

RETURN 

DEBUG  SUBCHK 
END 
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•OPTIONS  IN  EFFECT*  NOTERM. I O.EBCOI C . 
•OPTIONS  IN  EFFECT*  NAME  ■ MATOUT  , 
•STATISTICS*  SOURCE  STATEMENTS  * 
•STATISTICS*  NO  DIAGNOSTICS  GENERATED 


SOURCE • NOL I ST • NODECK. LOAD. NOMAP. NOTEST 
LINECNT  m 60 

11 .PROGRAM  SIZE  • 700 


•STATISTICS* 


NO  DIAGNOSTICS  THIS  STEP 
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INITIAL  DISTRIBUTION 


HQ  USAF/SAMI  1 

USAFE/DOQ  1 

PACAF/DOOFQ  1 

TAC/DRA  1 

ASD/ENFEA  1 

AUL/LSE-  71-249  1 

SAC/NRI  (STINFO  LIB)  1 

NWC/CODE  318  1 

NWC/CODE  317  1 

OO-ALC/NMWMP  2 

AFIS/INTA  1 

DDC  2 

AFATL/DLODL  2 

AFATL/DL  1 

AFATL/DLY  1 

ADTC/XRS  1 

AFATL/DLYV  20 

AFATL/DLYW  10 

USA  ENG  WatWay  Ex  Sta/VMS  1 

BAL  RESRCH  LAB/AMXBR-VL  1 

AMSAA/DRXSY-J  2 

USAMSAA/DRXSY-S  1 

ARRADCCM/DRDAR-LCU-1M  1 

AFOSR/N4  1 

ARMY  RESEARCH  OFFICE/NC  1 

OKLAHOMA  ST  UNIV/Dept  of  Stat  20 
TAC/INAT  1 

USA  TRADOC  SYS  ANN  ACT  1 

ASD/XRP  1 

COMIPAC/I-232  1 
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